\section{随机模拟积分}

某些非随机的问题也可以通过概率模型引入随机变量, 化为求随机模型的未知参数的问题。§3.1中用随机投点法估计 \(\pi\) 就是这样的一个例子。

{\color{red}\textbf{[核心思想]:} 随机模拟积分的本质是``用随机性解决确定性问题''——把计算积分转化为求期望,把几何问题转化为概率问题。这是一种反直觉但极其强大的思想。}

随机模拟解决非随机问题的典型代表是计算定积分。通过对随机模拟定积分的讨论, 可以展示随机模拟中大部分的问题和技巧。随机模拟积分也称为蒙特卡洛 (Monte Carlo) 积分 (简称 MC 积分)。实际上, 统计中最常见的计算问题就是积分和最优化问题。

\subsection{随机投点法}

设函数 \(h\left( x\right)\) 在有限区间 \(\left\lbrack  {a,b}\right\rbrack\) 上定义且有界,不妨设 \(0 \leq  h\left( x\right)  \leq  M\) 。要计算 \(I =\)  \({\int }_{a}^{b}h\left( x\right) {dx}\) ,相当于计算曲线下的区域 \(D = \{ \left( {x,y}\right)  : 0 \leq  y \leq  h\left( x\right) ,x \in  C = \left\lbrack  {a,b}\right\rbrack  \}\) 的面积。

{\color{red}\textbf{[物理图像]:} 想象向一个矩形靶子随机投飞镖,曲线下方区域占总靶面积的比例,就等于飞镖落入该区域的概率。投得越多,频率越接近真实比例,面积也就算得越准。}

为此在 \(G = \left\lbrack  {a,b}\right\rbrack   \times  \left( {0,M}\right)\) 上均匀抽样 \(N\) 次,得随机点 \({Z}_{1},{Z}_{2}\ldots ,{Z}_{N},{Z}_{i} = \left( {{X}_{i},{Y}_{i}}\right) ,i = 1,2,\ldots ,N\) 。 令

\[
{\xi }_{i} = \left\{  {\begin{array}{ll} 1, & {Z}_{i} \in  D \\  0, & \text{ 其它 } \end{array},\;i = 1,\ldots ,N}\right.
\]

则 \(\left\{  {\xi }_{i}\right\}\) 是独立重复试验结果, \(\left\{  {\xi }_{i}\right\}\) 独立同 \(\mathrm{b}\left( {1,p}\right)\) 分布,

\[
p = P\left( {{Z}_{i} \in  D}\right)  = V\left( D\right) /V\left( G\right)  = I/\left\lbrack  {M\left( {b - a}\right) }\right\rbrack   \tag{3.2}
\]

其中 \(V\left( \cdot \right)\) 表示区域面积。

\begin{center}
\includegraphics[max width=0.8\textwidth]{images/bo_d3dpgvk601uc738iknq0_3_311_866_1126_762_0.jpg}
\end{center}
\hspace*{3em} 

图 3.1: 随机投点法积分示意图

从模拟产生的随机样本 \({Z}_{1},{Z}_{2}\ldots ,{Z}_{N}\) ,可以用这 \(N\) 个点中落入曲线下方区域 \(D\) 的百分比 \(\widehat{p}\) 来估计 (3.2)中的概率 \(p\) (见图3.1),然后由 \(I = {pM}\left( {b - a}\right)\) 得到定积分 \(I\) 的近似值

\[
\widehat{I} = \widehat{p}M\left( {b - a}\right)  \tag{3.3}
\]

这种方法叫做随机投点法。这样计算的定积分有随机性, 误差中包含了随机模拟误差。

由强大数律可知

\[
\widehat{p} = \frac{\sum {\xi }_{i}}{N} \rightarrow  p,\text{ a.s. }\left( {N \rightarrow  \infty }\right)
\]

\[
\widehat{I} = \widehat{p}M\left( {b - a}\right)  \rightarrow  {pM}\left( {b - a}\right)  = I\text{ , a.s. }\left( {N \rightarrow  \infty }\right)
\]

即 \(N \rightarrow  \infty\) 时精度可以无限地提高 (当然,在计算机中要受到数值精度的限制)。

那么, 提高精度需要多大的代价呢? 由中心极限定理可知

\[
\sqrt{N}\left( {\widehat{p} - p}\right) /\sqrt{p\left( {1 - p}\right) }\overset{\mathrm{d}}{ \rightarrow  }\mathrm{N}\left( {0,1}\right) ,\left( {N \rightarrow  \infty }\right) ,
\]

从而

\[
\sqrt{N}\left( {\widehat{I} - I}\right)  = M\left( {b - a}\right) \left( {\widehat{p} - p}\right) \overset{\mathrm{d}}{ \rightarrow  }\mathrm{N}\left( {0,{\left\lbrack  M\left( b - a\right) \right\rbrack  }^{2}p\left( {1 - p}\right) }\right)  \tag{3.4}
\]

{\color{red}\textbf{[数学洞察]:} 这个$O_p(1/\sqrt{N})$的收敛速度是随机模拟的``魔咒''——精度提高10倍需要样本量增加100倍!这是用随机性付出的代价,但在高维问题中这个代价反而成为优势。}

当 \(N\) 很大时 \(\widehat{I}\) 近似服从 \(\mathrm{N}\left( {I,{\left\lbrack  M\left( b - a\right) \right\rbrack  }^{2}p\left( {1 - p}\right) /N}\right)\) 分布,称此近似分布的方差 \(\lbrack M(b -\)  \(\left. {\left. a\right) {\rbrack }^{2}p\left( {1 - p}\right) /N}\right)\) 为 \(\widehat{I}\) 的渐近方差。计算渐近方差可以用 \(\widehat{p}\) 代替 \(p\) 估计为 \({\left\lbrack  M\left( b - a\right) \right\rbrack  }^{2}\widehat{p}\left( {1 - \widehat{p}}\right) /N\) 。 (3.4)说明 \(\widehat{I}\) 的误差为 \({O}_{p}\left( \frac{1}{\sqrt{N}}\right)\) ,这样,计算 \(\widehat{I}\) 的精度每增加一位小数,计算量需要增加 100 倍。随机模拟积分一般都服从这样的规律。

\subsection{平均值法}

为了计算 \(I = {\int }_{a}^{b}h\left( x\right) {dx}\) ,上面用了类似于 \(§{2.2.4}\) 的舍选法的做法,在非随机问题中引入随机性时用了二维均匀分布和两点分布,靠求两点分布概率来估计积分 \(I\) 。随机投点法容易理解,但是效率较低,另一种效率更高的方法是利用期望值的估计。取 \(U \sim  \mathrm{U}\left( {a,b}\right)\) ,则

\[
E\left\lbrack  {h\left( U\right) }\right\rbrack   = {\int }_{a}^{b}h\left( u\right) \frac{1}{b - a}{du} = \frac{I}{b - a}
\]

\[
I = \left( {b - a}\right)  \cdot  {Eh}\left( U\right)
\]

{\color{red}\textbf{[关键转变]:} 从``投点计数''到``函数求平均''——平均值法不需要额外的$Y$坐标,只在$x$轴上随机取点,直接对函数值求平均。这减少了一个维度的随机性,效率自然更高。}

若取 \(\left\{  {{U}_{i},i = 1,\ldots ,N}\right\}\) 独立同 \(\mathrm{U}\left( {a,b}\right)\) 分布,则 \({Y}_{i} = h\left( {U}_{i}\right) ,i = 1,2,\ldots ,N\) 是 iid 随机变量列, 由强大数律,

\[
\bar{Y} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}h\left( {U}_{i}\right)  \rightarrow  {Eh}\left( U\right)  = \frac{I}{b - a},\;\text{ a.s. }\left( {N \rightarrow  \infty }\right)
\]

于是

\[
\widetilde{I} = \frac{b - a}{N}\mathop{\sum }\limits_{{i = 1}}^{N}h\left( {U}_{i}\right)  \tag{3.5}
\]

是 \(I\) 的强相合估计。称这样计算定积分 \(I\) 的方法为平均值法。由中心极限定理有

\[
\sqrt{N}\left( {\widetilde{I} - I}\right) \overset{\mathrm{d}}{ \rightarrow  }\mathrm{N}\left( {0,{\left( b - a\right) }^{2}\operatorname{Var}\left( {h\left( U\right) }\right) }\right)
\]

其中

\[
\operatorname{Var}\left\lbrack  {h\left( U\right) }\right\rbrack   = {\int }_{a}^{b}{\left\lbrack  h\left( u\right)  - Eh\left( U\right) \right\rbrack  }^{2}\frac{1}{b - a}{du} \tag{3.6}
\]

仅与 \(h\) 有关,仍有 \(\widetilde{I} - I = {O}_{p}\left( \frac{1}{\sqrt{N}}\right)\) ,但是(3.5)的渐近方差小于(3.3)的渐近方差 (见 \(§{3.2.3}\) )。 \(\operatorname{Var}\left\lbrack  {h\left( U\right) }\right\rbrack\) 可以用模拟样本 \(\left\{  {{Y}_{i} = h\left( {U}_{i}\right) }\right\}\) 估计为

\[
\operatorname{Var}\left( {h\left( U\right) }\right)  \approx  \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{\left( {Y}_{i} - \bar{Y}\right) }^{2}. \tag{3.7}
\]

如果定积分区间是无穷区间,比如 \({\int }_{0}^{\infty }h\left( x\right) {dx}\) ,为了使用均匀分布随机数以及平均值法计算积分可以做积分变换,使积分区间变成有限区间。例如,作变换 \(t = 1/\left( {x + 1}\right)\) ,则

\[
{\int }_{0}^{\infty }h\left( x\right) {dx} = {\int }_{0}^{1}h\left( {\frac{1}{t} - 1}\right) \frac{1}{{t}^{2}}{dt}.
\]

从平均值法看出,定积分问题 \({\int }_{a}^{b}h\left( x\right) {dx}\) 等价于求 \({Eh}\left( U\right)\) ,其中 \(U \sim  \mathrm{U}\left( {a,b}\right)\) 。所以这一节讨论的方法也是用来求随机变量函数期望的随机模拟方法。对一般随机变量 \(X\) ,其取值范围不必局限于有限区间,为了求 \(X\) 的函数 \(h\left( X\right)\) 的期望 \(I = {Eh}\left( X\right)\) ,对 \(X\) 的随机数 \({X}_{i},i = 1,2,\ldots ,N\) ,令 \({Y}_{i} = h\left( {X}_{i}\right)\) ,也可以用平均值法 \(\widehat{I} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}h\left( {X}_{i}\right)\) 来估计 \({Eh}\left( X\right) ,\widehat{I}\) 是 \(\operatorname{Eh}\left( X\right)\) 的无偏估计和强相合估计,若 \(\operatorname{Var}\left( {h\left( X\right) }\right)\) 存在,则 \(\widehat{I}\) 的方差为 \(\frac{1}{N}\operatorname{Var}\left( {h\left( X\right) }\right) ,\widehat{I}\) 有渐近正态分布 \(\mathrm{N}\left( {\operatorname{Eh}\left( X\right) ,\frac{1}{N}\operatorname{Var}\left( {h\left( X\right) }\right) }\right)\) 。设 \(\left\{  {Y}_{i}\right\}\) 的样本方差为 \({S}_{N}^{2}\) ,可以用 \({S}_{N}^{2}/N\) 估计 \(\widehat{I}\) 的方差,用 \(\widehat{I} \pm  2{S}_{N}/\sqrt{N}\) 作为 \(I\) 的近似 \({95}\%\) 置信区间。

例 3.2.1. 设 \(X \sim  \mathrm{N}\left( {0,1}\right)\) ,求 \(I = E{\left| X\right| }^{\frac{3}{2}}\) 。

解: 作变量替换积分可得

\[
I = 2{\int }_{0}^{\infty }{x}^{\frac{3}{2}}\frac{1}{\sqrt{2\pi }}{e}^{-\frac{1}{2}{x}^{2}}{dx}\;\left( {\text{ 令 }t = \frac{1}{2}{x}^{2}}\right)
\]

\[
= \frac{{2}^{\frac{3}{4}}\Gamma \left( \frac{5}{4}\right) }{\sqrt{\pi }} \approx  {0.86004}
\]

如果用平均值法估计 \(I\) ,取抽样样本量 \(N = {10000}\) ,产生标准正态随机数 \({X}_{i},i =\)  \(1,2,\ldots ,n\) ,令 \({Y}_{i} = {\left| {X}_{i}\right| }^{\frac{3}{2}}\) ,令 \(\widehat{I} = \bar{Y} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{\left| {X}_{i}\right| }^{\frac{3}{2}},{S}_{N}^{2} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{\left( {\left| {X}_{i}\right| }^{\frac{3}{2}} - \widehat{I}\right) }^{2}\) ,一次模拟的结果得到 \(\widehat{I} = {0.8658},{S}_{N} = {0.9409}\) ,渐近标准差为 \({S}_{N}/\sqrt{N} = {0.00941}\) ,说明 \({95}\%\) 置信水平下误差界为 \(\pm  {0.02}\) ,结果只有将近两位小数的精度。为了达到四位小数的精度,需要误差控制在 \(\pm  {0.00005}\) 以下,需要 \(2{S}_{N}/\sqrt{N} < {0.00005}\) ,代入得 \(N\) 需要超过 90 亿,可见随机模拟方法提高精度的困难程度。

以上的评估误差的方法使用了渐近方差和渐近正态分布, 在有些模拟问题中估计量不一定有渐近方差, 或渐近分布不一定准确 (比如产生的随机数不独立时)。在计算量允许的情形下,可以用不同的随机数种子重复得到 \(B\) 个 \(\widehat{I}\) 的估计值,用这 \(B\) 个 \(\widehat{I}\) 的样本标准差来估计 \(\widehat{I}\) 的抽样分布标准差。取 \(B = {100}\) 的一次计算结果得到的 \(\widehat{I}\) 的抽样分布标准差为 0.00924, 与从一次模拟得到的渐近标准差结果十分接近。也可以使用 \(§{3.6}\) 的 bootstrap 方法获取 \(B\) 组 bootstrap 样本,得到 \(B\) 个 bootstrap 的 \(\widehat{I}\) 样本值,从中估计抽样分布标准差。bootstrap 方法可以避免重新计算 \(h\left( {X}_{i}\right)\) 的值,在更复杂的问题中往往 \(h\left( {X}_{i}\right)\) 的计算是速度很慢的,比如说, \(h\left( {X}_{i}\right)\) 本身也需要用随机模拟或者数值优化方法计算。

实际上,一维积分用数值方法均匀布点计算一般更有效。比如,令 \({x}_{i} = a + \frac{b - a}{N}i\) , \(i = 0,1,\ldots ,N\) ,估计 \(I\) 为 (复合梯形求积公式)

\[
{I}_{0} = \frac{b - a}{N}\left\lbrack  {\frac{1}{2}h\left( a\right)  + \mathop{\sum }\limits_{{i = 1}}^{{N - 1}}h\left( {x}_{i}\right)  + \frac{1}{2}h\left( b\right) }\right\rbrack   \tag{3.8}
\]

{\color{red}\textbf{[对比]:} 确定性网格法精度$O(N^{-2})$远胜随机模拟的$O_p(N^{-1/2})$——对于光滑的一维积分,蒙特卡洛完全是``杀鸡用牛刀''。但到了高维,网格法指数爆炸,随机模拟才显神威。}

则当 \(h\left( x\right)\) 在 \(\left\lbrack  {a,b}\right\rbrack\) 上二阶连续可微时误差 \({I}_{0} - I\) 仅为 \(O\left( \frac{1}{{N}^{2}}\right)\) 阶,比随机模拟方法得到的精度要高得多,而且当 \(h\left( x\right)\) 有更好的光滑性时还可以用更精确的求积公式得到更高精度。所以, 被积函数比较光滑的一元定积分问题一般不需要用随机模拟来计算。

\subsection{高维定积分}

上面的两种计算一元函数定积分的方法可以很容易地推广到多元函数定积分, 或称高维定积分。设 \(d\) 元函数 \(h\left( {{x}_{1},{x}_{2},\ldots ,{x}_{d}}\right)\) 定义于超矩形

\[
C = \left\{  {\left( {{x}_{1},{x}_{2},\ldots ,{x}_{d}}\right)  : {a}_{i} \leq  {x}_{i} \leq  {b}_{i},i = 1,2,\ldots ,d}\right\}
\]

且

\[
0 \leq  h\left( {{x}_{1},\ldots ,{x}_{d}}\right)  \leq  M,\forall x \in  C.
\]

令

\[
D = \left\{  {\left( {{x}_{1},{x}_{2},\ldots ,{x}_{d},y}\right)  : \left( {{x}_{1},{x}_{2},\ldots ,{x}_{d}}\right)  \in  C,0 \leq  y \leq  h\left( {{x}_{1},{x}_{2},\ldots ,{x}_{d}}\right) }\right\}  ,
\]

\[
G = \left\{  {\left( {{x}_{1},{x}_{2},\ldots ,{x}_{d},y}\right)  : \left( {{x}_{1},{x}_{2},\ldots ,{x}_{d}}\right)  \in  C,0 \leq  y \leq  M}\right\}
\]

为计算 \(d\) 维定积分

\[
I = {\int }_{{a}_{d}}^{{b}_{d}}\cdots {\int }_{{a}_{2}}^{{b}_{2}}{\int }_{{a}_{1}}^{{b}_{1}}h\left( {{x}_{1},{x}_{2},\ldots ,{x}_{d}}\right) d{x}_{1}d{x}_{2}\cdots d{x}_{d}, \tag{3.9}
\]

产生服从 \(d + 1\) 维空间中的超矩形 \(G\) 内的均匀分布的独立抽样 \({\mathbf{Z}}_{1},{\mathbf{Z}}_{2},\ldots ,{\mathbf{Z}}_{N}\) ,令

\[
{\xi }_{i} = \left\{  {\begin{array}{ll} 1, & {\mathbf{Z}}_{i} \in  D \\  0, & {\mathbf{Z}}_{i} \in  G - D \end{array},\;i = 1,2,\ldots ,N}\right.
\]

则 \({\xi }_{i}\) iid \(\mathrm{b}\left( {1,p}\right)\) ,

\[
p = P\left( {{\mathbf{Z}}_{i} \in  D}\right)  = \frac{V\left( D\right) }{V\left( G\right) } = \frac{I}{{MV}\left( C\right) } = \frac{I}{M\mathop{\prod }\limits_{{j = 1}}^{d}\left( {{b}_{j} - {a}_{j}}\right) }
\]

其中 \(V\left( \cdot \right)\) 表示区域体积。令 \(\widehat{p}\) 为 \(N\) 个随机点中落入 \(D\) 的百分比,则

\[
\widehat{p} = \frac{\sum {\xi }_{i}}{N} \rightarrow  p,\text{ a.s. }\left( {N \rightarrow  \infty }\right) ,
\]

用

\[
{\widehat{I}}_{1} = \widehat{p}V\left( G\right)  = \widehat{p} \cdot  {MV}\left( C\right)  = \widehat{p} \cdot  M\mathop{\prod }\limits_{{j = 1}}^{d}\left( {{b}_{j} - {a}_{j}}\right)  \tag{3.10}
\]

估计积分 \(I\) ,则 \({\widehat{I}}_{1}\) 是 \(I\) 的无偏估计和强相合估计。称用(3.10)计算高维定积分 \(I\) 的方法为随机投点法。由中心极限定理知

\[
\sqrt{N}\left( {\widehat{p} - p}\right) /\sqrt{p\left( {1 - p}\right) }\overset{d}{ \rightarrow  }\mathrm{\;N}\left( {0,1}\right)
\]

\[
\sqrt{N}\left( {{\widehat{I}}_{1} - I}\right) \overset{d}{ \rightarrow  }\mathrm{\;N}\left( {0,{\left( M\mathop{\prod }\limits_{{j = 1}}^{n}\left( {b}_{j} - {a}_{j}\right) \right) }^{2}p\left( {1 - p}\right) }\right) ,
\]

\({\widehat{I}}_{1}\) 的渐近方差为

\[
\frac{{\left( M\mathop{\prod }\limits_{{j = 1}}^{n}\left( {b}_{j} - {a}_{j}\right) \right) }^{2}p\left( {1 - p}\right) }{N} \tag{3.11}
\]

所以 \({\widehat{I}}_{1}\) 的随机误差仍为 \({O}_{p}\left( \frac{1}{\sqrt{N}}\right) ,N \rightarrow  \infty\) 时的误差阶不受维数 \(d\) 的影响,这是随机模拟方法与其它数值计算方法相比一个重大优势。

{\color{red}\textbf{[为什么高维有优势?]:} 网格法在$d$维需要$n^d$个点才能达到$O(n^{-2})$精度,即$O(N^{-2/d})$。当$d=10$时精度仅$O(N^{-0.2})$,而蒙特卡洛始终保持$O(N^{-0.5})$!这就是打破``维数诅咒''的关键。}

在计算高维积分(3.9)时,仍可以通过估计 \({Eh}\left( \mathbf{U}\right)\) 来获得,其中 \(\mathbf{U}\) 服从 \({R}^{d}\) 中的超矩形 \(C\) 上的均匀分布。设 \({\mathbf{U}}_{i} \sim\) iid \(\mathrm{U}\left( C\right) ,i = 1,2,\ldots ,N\) ,则 \(h\left( {\mathbf{U}}_{i}\right) ,i = 1,2,\ldots ,N\) 是 iid 随机变量列,

\[
{Eh}\left( {\mathbf{U}}_{i}\right)  = {\int }_{C}h\left( \mathbf{u}\right) \frac{1}{\mathop{\prod }\limits_{{j = 1}}^{d}\left( {{b}_{j} - {a}_{j}}\right) }d\mathbf{u} = \frac{I}{\mathop{\prod }\limits_{{j = 1}}^{d}\left( {{b}_{j} - {a}_{j}}\right) },
\]

估计 \(I\) 为

\[
{\widehat{I}}_{2} = \mathop{\prod }\limits_{{j = 1}}^{d}\left( {{b}_{j} - {a}_{j}}\right)  \cdot  \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}h\left( {U}_{i}\right) , \tag{3.12}
\]

称用(3.12)计算高维定积分 \(I\) 的方法为平均值法。由强大数律

\[
{\widehat{I}}_{2} \rightarrow  \mathop{\prod }\limits_{{j = 1}}^{d}\left( {{b}_{j} - {a}_{j}}\right) {Eh}\left( U\right)  = I,\text{ a.s. }\left( {N \rightarrow  \infty }\right) ,
\]

\({\widehat{I}}_{2}\) 的渐近方差为

\[
\frac{{\left( V\left( C\right) \right) }^{2}\operatorname{Var}\left( {h\left( \mathbf{U}\right) }\right) }{N} = \frac{{\left( \mathop{\prod }\limits_{{j = 1}}^{d}\left( {b}_{j} - {a}_{j}\right) \right) }^{2}\operatorname{Var}\left( {h\left( \mathbf{U}\right) }\right) }{N}. \tag{3.13}
\]

\(N \rightarrow  \infty\) 时的误差阶也不受维数 \(d\) 的影响。

我们来比较随机投点法(3.10)与平均值法(3.12)的精度,只要比较其渐近方差。对 \(I =\)  \({\int }_{C}h\left( \mathbf{x}\right) d\mathbf{x}\) ,设 \({\widehat{I}}_{1}\) 为随机投点法的估计, \({\widehat{I}}_{2}\) 为平均值法的估计。因设 \(0 \leq  h\left( \mathbf{x}\right)  \leq  M\) ,不妨设 \(0 \leq  h\left( \mathbf{x}\right)  \leq  1\) ,取 \(h\left( \mathbf{x}\right)\) 的上界 \(M = 1\) 。

令 \({\mathbf{X}}_{i} \sim\) iid \(\mathrm{U}\left( C\right) ,{\eta }_{i} = h\left( {\mathbf{X}}_{i}\right) ,{Y}_{i} \sim\) iid \(\mathrm{U}\left( {0,1}\right)\) 与 \(\left\{  {\mathbf{X}}_{i}\right\}\) 独立,

\[
{\xi }_{i} = \left\{  {\begin{array}{ll} 1, & \text{ 当 }{Y}_{i} \leq  h\left( {\mathbf{X}}_{i}\right) \\  0, & \text{ 当 }{Y}_{i} > h\left( {\mathbf{X}}_{i}\right)  \end{array}i = 1,2,\ldots ,N,}\right.
\]

这时有

\[
{\widehat{I}}_{1} = V\left( C\right) \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{\xi }_{i},\;{\widehat{I}}_{2} = V\left( C\right) \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{\eta }_{i}
\]

\[
\operatorname{Var}\left( {\widehat{I}}_{1}\right)  = \frac{1}{N}{V}^{2}\left( C\right)  \cdot  \frac{I}{V\left( C\right) }\left( {1 - \frac{I}{V\left( C\right) }}\right)
\]

\[
\operatorname{Var}\left( {\widehat{I}}_{2}\right)  = \frac{1}{N}{V}^{2}\left( C\right)  \cdot  \left( {\frac{1}{V\left( C\right) }{\int }_{C}{h}^{2}\left( \mathbf{x}\right) d\mathbf{x} - {\left( \frac{I}{V\left( C\right) }\right) }^{2}}\right)
\]

\[
\operatorname{Var}\left( {\widehat{I}}_{1}\right)  - \operatorname{Var}\left( {\widehat{I}}_{2}\right)  = \frac{V\left( C\right) }{N}{\int }_{C}\left\{  {h\left( \mathbf{x}\right)  - {h}^{2}\left( \mathbf{x}\right) }\right\}  d\mathbf{x} \geq  0
\]

可见平均值法精度更高。事实上,随机投点法多用了随机数 \({Y}_{i}\) ,必然会增加抽样随机误差。

{\color{red}\textbf{[解释]:} 投点法浪费在哪?它需要两个随机数$(X_i, Y_i)$判断是否在曲线下,而平均值法只需一个$X_i$直接计算函数值。多一个随机维度就多一份方差,这是效率差异的根源。}

在计算高维积分(3.10)时,如果用网格法作数值积分,把超矩形 \(C = \left\lbrack  {{a}_{1},{b}_{1}}\right\rbrack   \times  \left\lbrack  {{a}_{2},{b}_{2}}\right\rbrack   \times\)  \(\cdots  \times  \left\lbrack  {{a}_{d},{b}_{d}}\right\rbrack\) 的每个边分成 \(n\) 个格子,就需要 \(N = {n}^{d}\) 个格子点,如果用每个小超矩形的中心作为代表点,可以达到 \(O\left( {n}^{-2}\right)\) 精度,即 \(O\left( {N}^{-2/d}\right)\) ,当维数增加时为了提高一倍精度需要 \({2}^{d/2}\) 倍的代表点。比如 \(d = 8\) ,精度只有 \(O\left( {N}^{-1/4}\right)\) 。高维的问题当维数增加时计算量会迅猛增加,以至于无法计算,这个问题称为维数诅咒。

{\color{red}\textbf{[量级估计]:} 对$d=8$的问题,网格法要达到$10^{-4}$精度需要$(10^4)^{4}=10^{16}$个点!而蒙特卡洛只需$(10^4)^2=10^8$个点。这是$10^8$倍的差距,网格法根本无法计算。}

如果用 Monte Carlo 积分,则精度为 \({O}_{p}\left( {N}^{-1/2}\right)\) ,与 \(d\) 关系不大。所以 Monte Carlo 方法在高维积分中有重要应用。为了提高积分计算精度,需要减小 \({O}_{\mathrm{p}}\left( {N}^{-1/2}\right)\) 中的常数项,即减小 \(\widehat{I}\) 的渐近方差。

例 3.2.2. 考虑

\[
f\left( {{x}_{1},{x}_{2},\ldots ,{x}_{d}}\right)  = {x}_{1}^{2}{x}_{2}^{2}\cdots {x}_{d}^{2},0 \leq  {x}_{1} \leq  1,0 \leq  {x}_{2} \leq  1,\ldots ,0 \leq  {x}_{d} \leq  1
\]

的积分

\[
I = {\int }_{0}^{1}\cdots {\int }_{0}^{1}{\int }_{0}^{1}{x}_{1}^{2}{x}_{2}^{2}\cdots {x}_{d}^{2}d{x}_{1}d{x}_{2}\cdots d{x}_{d}
\]

当然,这个积分是有精确解 \(I = {\left( 1/3\right) }^{d}\) 的,对估计 \(I\) 的结果我们可以直接计算误差。以 \(d = 8\) 为例比较以下三种方法的精度,这时真值 \(I = {\left( 1/3\right) }^{8} \approx  {1.524} \times  {10}^{-4}\) 。

用 \(n\) 网格点法, \(N = {n}^{d}\) ,公式为

\[
{I}_{0} = \frac{1}{N}\mathop{\sum }\limits_{{{i}_{1} = 1}}^{n}\mathop{\sum }\limits_{{{i}_{2} = 1}}^{n}\ldots \mathop{\sum }\limits_{{{i}_{d} = 1}}^{n}f\left( {\frac{2{i}_{1} - 1}{2n},\frac{2{i}_{2} - 1}{2n},\ldots ,\frac{2{i}_{d} - 1}{2n}}\right)
\]

误差绝对值为 \({e}_{0} = \left| {{I}_{0} - I}\right|\) 。如果取 \(n = 5,d = 8\) ,需要计算 \(N = {390625}\) 个点的函数值,计算量相当大。

用随机投点法求 \(I\) ,先在 \({\left( 0,1\right) }^{d} \times  \left( {0,1}\right)\) 均匀抽样 \(\left( {{\xi }_{i}^{\left( 1\right) },{\xi }_{i}^{\left( 2\right) },\ldots ,{\xi }_{i}^{\left( d\right) },{y}_{i}}\right) ,i = 1,\ldots ,N\) ,令 \({\widehat{I}}_{1}\) 为 \({y}_{i} \leq  f\left( {{\xi }_{i}^{\left( 1\right) },{\xi }_{i}^{\left( 2\right) },\ldots ,{\xi }_{i}^{\left( d\right) }}\right)\) 成立的百分比。因为 \({\widehat{I}}_{1}\) 是随机的,误差绝对值 \(\left| {{\widehat{I}}_{1} - I}\right|\) 也是随机的,所以我们重复试验 \(B\) 次,计算 \(B\) 次的误差绝对值的平均值 \({e}_{1}\) ,作为 \(E\left| {{\widehat{I}}_{1} - I}\right|\) 的估计值。取 \(B\) 多大合适呢? 因为计算量很大,先取了 \({B}_{1} = {10}\) ,用得到的 \({B}_{1}\) 个 \(\left| {{\widehat{I}}_{1} - I}\right|\) 样本标准差 \({s}_{1}\) 可以估计 \(B\) 个 \(\left| {{\widehat{I}}_{1} - I}\right|\) 的样本平均值的标准差为 \({\mathrm{{SE}}}_{1} = {s}_{1}/\sqrt{B}\) ,发现 \(B = {1000}\) 时可以控制在相对误差 2\% 以下。

用平均值方法, 公式为

\[
{\widehat{I}}_{2} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}f\left( {{\xi }_{i}^{\left( 1\right) },{\xi }_{i}^{\left( 2\right) },\cdots ,{\xi }_{i}^{\left( d\right) }}\right)
\]

其中 \({\xi }_{i}^{\left( j\right) },i = 1,\ldots ,N,j = 1,\ldots ,d\) 独立同 \(\mathrm{U}\left( {0,1}\right)\) 分布。类似地,重复 \(B = {1000}\) 次,计算 \(B\) 次的误差绝对值 \(\left| {{\widehat{I}}_{2} - I}\right|\) 的平均值 \({e}_{2}\) ,作为 \(E\left| {{\widehat{I}}_{2} - I}\right|\) 的估计。

最后,三种方法计算量相同(都计算 \(N = {5}^{8}\) 次函数值)的情况下,得到网格点法的误差绝对值为 \({e}_{0} = {1.2} \times  {10}^{-5}\) ,随机投点法的误差绝对值平均值为 \({e}_{1} = {1.6} \times  {10}^{-5}\) ,平均值法的误差绝对值平均值为 \({e}_{2} = {0.20} \times  {10}^{-5}\) ,此问题结果中平均值法比其它两种方法精度高了至少 5 倍。

\subsection{重要抽样法}

\(I = {\int }_{C}h\left( \mathbf{x}\right) d\mathbf{x}\) 中积分区域 \(C\) 可能是任意形状的,也可能无界, \(h\left( \mathbf{x}\right)\) 在 \(C\) 内各处的取值大小差异可能很大,使得直接用平均值法估计 \(I\) 时,很多样本点处于 \(\left| {h\left( \mathbf{x}\right) }\right|\) 接近于零的地方, 造成浪费,另外使得 \({\widehat{I}}_{2}\) 的渐近方差 (见(3.13)) 中的 \(\operatorname{Var}\left( {h\left( \mathbf{U}\right) }\right)\) 很大 \(\left( {\mathbf{U} \sim  \mathrm{U}\left( C\right) }\right)\) 。

{\color{red}\textbf{[类比]:} 重要抽样就像``有的放矢''——如果被积函数只在山峰处有值、平原处为零,均匀撒点就是浪费弹药。聪明的做法是在山峰密集投点、平原稀疏投点,然后用权重补偿这种不均匀性。}

为此,考虑用非均匀抽样: \(\left| {h\left( \mathbf{x}\right) }\right|\) 大的地方密集投点, \(\left| {h\left( \mathbf{x}\right) }\right|\) 小的地方稀疏投点。这样可以有效利用样本。设 \(g\left( \mathbf{x}\right) ,\mathbf{x} \in  C\) 是一个密度,其形状与 \(\left| {h\left( \mathbf{x}\right) }\right|\) 相近,且当 \(g\left( \mathbf{x}\right)  = 0\) 时 \(h\left( \mathbf{x}\right)  = 0\) ,当 \(\parallel \mathbf{x}\parallel  \rightarrow  \infty\) 时 \(h\left( \mathbf{x}\right)  = o\left( {g\left( \mathbf{x}\right) }\right)\) 。称 \(g\left( \mathbf{x}\right)\) 为试投密度或重要抽样密度。

设 \({\mathbf{X}}_{i}\) iid \(\sim  g\left( \mathbf{x}\right) ,i = 1,2,\ldots ,N\) 。令

\[
{\eta }_{i} = \frac{h\left( {\mathbf{X}}_{i}\right) }{g\left( {\mathbf{X}}_{i}\right) },i = 1,2,\ldots ,N
\]

则

\[
E{\eta }_{1} = {\int }_{C}\frac{h\left( \mathbf{x}\right) }{g\left( \mathbf{x}\right) }g\left( \mathbf{x}\right) d\mathbf{x} = {\int }_{C}h\left( \mathbf{x}\right) d\mathbf{x} = I
\]

{\color{red}\textbf{[数学洞察]:} 重要抽样的核心公式$E[h(X)/g(X)] = \int h$看似魔术,实则是测度变换的直接应用——从$g$的测度下计算$h/g$的期望,等价于从Lebesgue测度下计算$h$的积分。}

因此可以用 \(\left\{  {{\eta }_{i},i = 1,2,\ldots ,N}\right\}\) 的样本平均值来估计 \(I\) ,即

\[
{\widehat{I}}_{3} = \bar{\eta } = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}\frac{h\left( {\mathbf{X}}_{i}\right) }{g\left( {\mathbf{X}}_{i}\right) }. \tag{3.14}
\]

\({\widehat{I}}_{3}\) 是 \(I\) 的无偏估计和强相合估计。 \({\widehat{I}}_{3}\) 的渐近方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{3}\right)  = \operatorname{Var}\left( \frac{h\left( \mathbf{X}\right) }{g\left( \mathbf{X}\right) }\right) \frac{1}{N} = O\left( \frac{1}{N}\right) , \tag{3.15}
\]

{\color{red}\textbf{[方差推导]:} 这个等式基于蒙特卡洛估计的方差性质。设$\eta_i = h(X_i)/g(X_i)$为独立同分布的随机变量,则样本均值$\widehat{I}_3 = \frac{1}{N}\sum_{i=1}^N \eta_i$的方差为:
\[
\operatorname{Var}(\widehat{I}_3) = \operatorname{Var}\left(\frac{1}{N}\sum_{i=1}^N \eta_i\right) = \frac{1}{N^2} \cdot N \cdot \operatorname{Var}(\eta_i) = \frac{1}{N}\operatorname{Var}\left(\frac{h(\mathbf{X})}{g(\mathbf{X})}\right).
\]
由于$\operatorname{Var}(h(\mathbf{X})/g(\mathbf{X}))$是不依赖于$N$的常数(假设有限),因此$\operatorname{Var}(\widehat{I}_3) = O(1/N)$,表示方差以$1/N$的速率递减。这是蒙特卡洛方法的标准收敛速率。}

当 \(g\left( \mathbf{x}\right)\) 和 \(\left| {h\left( \mathbf{x}\right) }\right|\) 形状很接近时 \(\frac{\left| h\left( \mathbf{x}\right) \right| }{g\left( \mathbf{x}\right) }\) 近似为常数,方差 \(\operatorname{Var}\left( {\widehat{I}}_{3}\right)\) 很小,这时 \({\widehat{I}}_{3}\) 的随机误差可以很小。

{\color{red}\textbf{[理想情况]:} 最优试投密度是$g(\mathbf{x}) \propto |h(\mathbf{x})|$,此时$h/g$确实为常数,方差为零!当然这需要知道$\int |h|$,相当于已经解决了问题。实际中我们只能用近似的$g$,但仍能大幅降低方差。}

用(3.14)估计 \(I\) 与 \(§{2.2.4}\) 的舍选法 II 有类似的想法,这种方法叫做重要抽样法 (importance sampling), 是随机模拟的重要方法。

为什么当 \(g\left( \mathbf{x}\right)\) 和 \(\left| {h\left( \mathbf{x}\right) }\right|\) 形状很接近时 \(\operatorname{Var}\left( {\widehat{I}}_{3}\right)\) 很小？事实上,

{\color{red}\textbf{[证明目标]:} 我们要严格证明:方差$\operatorname{Var}(\widehat{I}_3)$有下界,且当$g(\mathbf{x}) \propto |h(\mathbf{x})|$时达到最小值。证明策略是利用Jensen不等式建立方差的下界。}

\[
\operatorname{Var}\left( \frac{h\left( \mathbf{X}\right) }{g\left( \mathbf{X}\right) }\right)  = E\left( \frac{{h}^{2}\left( \mathbf{X}\right) }{{g}^{2}\left( \mathbf{X}\right) }\right)  - \left( {E\frac{h\left( \mathbf{X}\right) }{g\left( \mathbf{X}\right) }}\right)  = E\left( \frac{{h}^{2}\left( \mathbf{X}\right) }{{g}^{2}\left( \mathbf{X}\right) }\right)  - {\left( {\int }_{C}h\left( \mathbf{x}\right) d\mathbf{x}\right) }^{2}
\]

{\color{red}\textbf{[为什么用Jensen不等式]:} 注意到$E(h^2/g^2)$很难直接处理,但我们可以利用$h^2 = (|h|)^2$,然后对凸函数$\phi(t) = t^2$应用Jensen不等式。Jensen不等式说:若$\phi$是凸函数,则$E[\phi(Y)] \geq \phi(E[Y])$。这里取$Y = |h(\mathbf{X})|/g(\mathbf{X})$,则$\phi(Y) = Y^2 = h^2/g^2$,于是$E(h^2/g^2) \geq (E|h|/g)^2$。}

{\color{red}\textbf{[数学推导第一步]:} 应用Jensen不等式:由于$t \mapsto t^2$是凸函数,对随机变量$|h(\mathbf{X})|/g(\mathbf{X})$有
\[E\left(\frac{h^2(\mathbf{X})}{g^2(\mathbf{X})}\right) = E\left[\left(\frac{|h(\mathbf{X})|}{g(\mathbf{X})}\right)^2\right] \geq \left[E\left(\frac{|h(\mathbf{X})|}{g(\mathbf{X})}\right)\right]^2.\]
因此方差$\geq (E|h|/g)^2 - I^2$。}

\[
\geq  {\left( E\frac{\left| h\left( \mathbf{X}\right) \right| }{g\left( \mathbf{X}\right) }\right) }^{2} - {I}^{2}\;\text{ (由 Jensen 不等式) } \tag{3.16}
\]

{\color{red}\textbf{[数学推导第二步]:} 计算$E(|h(\mathbf{X})|/g(\mathbf{X}))$:注意$\mathbf{X} \sim g(\mathbf{x})$,所以
\[E\left(\frac{|h(\mathbf{X})|}{g(\mathbf{X})}\right) = \int_C \frac{|h(\mathbf{x})|}{g(\mathbf{x})} \cdot g(\mathbf{x}) d\mathbf{x} = \int_C |h(\mathbf{x})| d\mathbf{x}.\]
这就是被积函数绝对值的积分!}

\[
= {\left( {\int }_{C}\left| h\left( \mathbf{x}\right) \right| d\mathbf{x}\right) }^{2} - {I}^{2}, \tag{3.17}
\]

{\color{red}\textbf{[最优条件分析]:} Jensen不等式取等号的条件是随机变量$|h(\mathbf{X})|/g(\mathbf{X})$为常数(几乎处处)。设这个常数为$k$,即$|h(\mathbf{x})| = k \cdot g(\mathbf{x})$对所有$\mathbf{x} \in C$成立。对两边积分:$\int_C |h(\mathbf{x})| d\mathbf{x} = k \int_C g(\mathbf{x}) d\mathbf{x} = k$。因此$k = \int_C |h(\mathbf{x})| d\mathbf{x}$,所以最优试投密度为$g(\mathbf{x}) = |h(\mathbf{x})|/k = |h(\mathbf{x})|/\int_C |h(\mathbf{t})| dt$。}

{\color{red}\textbf{[直觉理解]:} 为什么$g \propto |h|$最优?因为此时每个样本点$\mathbf{X}_i$都``等价重要'':虽然$h(\mathbf{X}_i)/g(\mathbf{X}_i)$的绝对值是常数,但符号可能不同(如果$h$有正有负)。方差来自$h$的符号变化,而不是幅度变化。如果$h \geq 0$,则$h/g$完全为常数,方差为零!}

当且仅当 \(\frac{\left| h\left( \mathbf{X}\right) \right| }{g\left( \mathbf{X}\right) }\) 为常数时不等式(3.16)中的等号成立,即 \(g\left( \mathbf{x}\right)  = \frac{1}{{\int }_{C}\left| {h\left( \mathbf{t}\right) }\right| {dt}}\left| {h\left( \mathbf{x}\right) }\right| ,\mathbf{x} \in  C\) 时 \(\operatorname{Var}\left( {\widehat{I}}_{3}\right)\) 达到最小。

上述求积分的问题也可以表述为求随机变量函数的期望问题。设 \(\mathbf{Y} \sim  f\left( \mathbf{y}\right)\) ,要求 \(\mathbf{Y}\) 的函数 \(h\left( \mathbf{Y}\right)\) 的期望值 \({Eh}\left( \mathbf{Y}\right)  = \int h\left( \mathbf{y}\right) f\left( \mathbf{y}\right) d\mathbf{y}\) ,可以抽取 \(\mathbf{Y}\) 的随机数 \({\mathbf{Y}}_{i},i = 1,2,\ldots ,N\) ,然后用平均值法 \(\frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}h\left( {\mathbf{Y}}_{i}\right)\) 估计 \({Eh}\left( \mathbf{Y}\right)\) 。但是,如果很难生成 \(\mathbf{Y}\) 的随机数,或者 \(\mathbf{Y}\) 的分布集中于 \(h\left( \mathbf{y}\right)\) 接近于零的位置以至于积分效率很低,可以找一个试投密度 \(g\left( \mathbf{x}\right)\) ,从 \(g\left( \mathbf{x}\right)\) 产生随机数 \({\mathbf{X}}_{i},i = 1,2,\ldots ,N\) ,用如下重要抽样法

\[
{\widehat{I}}_{3.1} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}h\left( {\mathbf{X}}_{i}\right) \frac{f\left( {\mathbf{X}}_{i}\right) }{g\left( {\mathbf{X}}_{i}\right) } \tag{3.18}
\]

估计 \({Eh}\left( \mathbf{Y}\right)\) ,易见 \(E{\widehat{I}}_{3.1} = {Eh}\left( \mathbf{Y}\right)\) 。称 \({W}_{i} = \frac{f\left( {\mathbf{X}}_{i}\right) }{g\left( {\mathbf{X}}_{i}\right) }\) 为重要性权重, \({Eh}\left( \mathbf{Y}\right)\) 可以用重要性权重估计为

\[
{\widehat{I}}_{3.1} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{W}_{i}h\left( {\mathbf{X}}_{i}\right) . \tag{3.19}
\]

选取试投密度 \(g\left( \mathbf{x}\right)\) 时,要求当 \(h\left( \mathbf{x}\right) f\left( \mathbf{x}\right)  \neq  0\) 时 \(g\left( \mathbf{x}\right)  \neq  0\) ,当 \(\mathbf{x}\) 趋于无穷时 \(h\left( \mathbf{x}\right) f\left( \mathbf{x}\right)  =\)  \(o\left( {g\left( \mathbf{x}\right) }\right)\) ( \(g\left( \mathbf{x}\right)\) 相对厚尾),一般还要求 \(g\left( \mathbf{x}\right)\) 的形状与 \(\left| {h\left( \mathbf{x}\right) }\right| f\left( \mathbf{x}\right)\) 形状接近。如果 \(g\left( \mathbf{x}\right)\) 的相对厚尾性难以确定, 可以使用如下保险的试投密度

\[
\widetilde{g}\left( \mathbf{x}\right)  = {\rho g}\left( \mathbf{x}\right)  + \left( {1 - \rho }\right) r\left( \mathbf{x}\right) , \tag{3.20}
\]

其中 \(\rho\) 接近于 1, \(r\left( \mathbf{x}\right)\) 是柯西分布或 Pareto 分布这样的重尾分布。要生成 \(N\) 个 \(\widetilde{g}\left( \mathbf{x}\right)\) 的随机数,只要生成 \({N\rho }\) 个 \(g\left( \mathbf{x}\right)\) 的随机数和 \(N\left( {1 - \rho }\right)\) 个 \(r\left( \mathbf{x}\right)\) 的随机数。

选取不适当的试投密度会把绝大多数样本点投到了对计算积分不重要的位置, 使得样本点中只有极少数点是真正有作用的。在多维问题中合适的试投密度尤其难找, 经常需要反复试验。

有时 \(\mathbf{Y} \sim  f\left( \cdot \right)\) 不仅很难直接抽样,而且 \(f\left( \cdot \right)\) 本身未知,只能确定到差一个常数倍的 \(\widetilde{f}\left( \mathbf{x}\right)  = {cf}\left( \mathbf{x}\right)\) ,常数 \(c\) 未知,为了求常数 \(c\) 需要计算 \(c = \int \widetilde{f}\left( \mathbf{y}\right) d\mathbf{y}\) ,计算 \(c\) 一般很困难。这时,定义重要性权重为 \({W}_{i} = \frac{\widetilde{f}\left( {\mathbf{X}}_{i}\right) }{g\left( {\mathbf{X}}_{i}\right) }\) ,公式(3.19)可以改成

\[
{\widehat{I}}_{4} = \frac{\mathop{\sum }\limits_{{i = 1}}^{N}{W}_{i}h\left( {\mathbf{X}}_{i}\right) }{\mathop{\sum }\limits_{{i = 1}}^{N}{W}_{i}} \tag{3.21}
\]

这称为标准化重要抽样法。对(3.21)的分子和分母都除以 \({cN}\) 后分子 a.s. 收敛到 \({Eh}\left( \mathbf{Y}\right)\) ,分母 a.s. 收敛到 1,所以 (3.21) 是 \(\operatorname{Eh}\left( \mathbf{Y}\right)\) 的强相合估计,但不是无偏的。标准化重要抽样估计往往比无偏估计 \({\widehat{I}}_{3.1}\) 有更小的均方误差。关于标准化重要抽样法的渐近方差的讨论参见 \(\operatorname{Liu}{\left( {2001}\right) }^{\left\lbrack  {28}\right\rbrack  }§{2.5.3}\) 。

如果需要对多个不同的函数 \(h\left( \cdot \right)\) 计算 \({Eh}\left( \mathbf{Y}\right)\) ,则选取试抽样密度 \(g\left( \mathbf{x}\right)\) 时应使得 \(g\left( \mathbf{x}\right)\) 尽可能与 \(\mathbf{Y}\) 的密度 \(f\left( \mathbf{x}\right)\) 形状接近,这样权重 \({W}_{i} = \frac{\widehat{f}\left( {\mathbf{X}}_{i}\right) }{g\left( {\mathbf{X}}_{i}\right) }\) 的分布不过于偏斜,不至于出现绝大部分权重集中于少数样本点的情形。抽样值 \({\mathbf{X}}_{i}\) 与权重 \({W}_{i}\) 一起可以看作是分布 \(f\left( \cdot \right)\) 的某种抽样。

{\color{red}\textbf{[背景说明]:} 前面的重要抽样中,样本点$\mathbf{X}_i$来自试投密度$g$,通过权重$W_i = f(\mathbf{X}_i)/g(\mathbf{X}_i)$来``纠正''分布。但实际中我们可能用更复杂的方法产生$({\mathbf{X}_i, W_i})$对(如MCMC、粒子滤波)。``适当加权抽样''定义了什么样的加权样本是``有效的''。}

定义 (适当加权抽样) 随机变量序列 \(\left\{  {\left( {{\mathbf{X}}_{i},{W}_{i}}\right) ,i = 1,2,\ldots ,N}\right\}\) 称为关于密度 \(f\left( \cdot \right)\) 的适当加权抽样,如果对于任何平方可积函数 \(h\left( \cdot \right)\) 都有

\[
E\left\lbrack  {h\left( {X}_{i}\right) {W}_{i}}\right\rbrack   = c{E}_{f}\left\lbrack  {h\left( X\right) }\right\rbrack   = c\int h\left( \mathbf{x}\right) f\left( \mathbf{x}\right) d\mathbf{x},i = 1,2,\ldots ,N,
\]

其中 \(c\) 是归一化常数。

{\color{red}\textbf{[定义的直觉理解]:} ``适当加权''的核心思想:虽然样本点$\mathbf{X}_i$不是直接从目标密度$f$抽取的,但通过权重$W_i$加权后,期望$E[h(\mathbf{X}_i)W_i]$与目标期望$E_f[h(X)]$成正比(相差常数$c$)。这意味着用$\sum h(\mathbf{X}_i)W_i / \sum W_i$可以一致地估计$E_f[h(X)]$,对所有函数$h$都成立!}

{\color{red}\textbf{[为什么需要常数$c$]:} 常数$c$的存在是因为:1) 实际中$f$可能只知道到差一个常数($\widetilde{f} = cf$),2) 权重$W_i$的尺度可以任意缩放而不影响标准化估计$\sum h(\mathbf{X}_i)W_i / \sum W_i$。关键是比例关系,而非绝对值。}

设随机变量(X, W)联合密度为 \(g\left( {\mathbf{x},w}\right)\) ,则(X, W)的样本为密度 \(f\left( \cdot \right)\) 的适当加权抽样的充分必要条件是

\[
{E}_{g}\left( {W \mid  \mathbf{x}}\right)  = {E}_{g}\left( W\right) \frac{f\left( \mathbf{x}\right) }{g\left( \mathbf{x}\right) },\forall \mathbf{x},
\]

其中 \({E}_{g}\left( W\right)\) 是关于 \(W\) 的边缘密度的期望, \({E}_{g}\left( {W \mid  \mathbf{x}}\right)\) 是在(X, W)的联合密度下条件期望 \(E\left( {W \mid  \mathbf{X}}\right)\) 在 \(\mathbf{X} = \mathbf{x}\) 处的值。

{\color{red}\textbf{[充要条件的含义]:} 这个条件说:给定位置$\mathbf{x}$,权重$W$的条件期望必须正比于似然比$f(\mathbf{x})/g(\mathbf{x})$,比例系数是$E_g(W)$。直观地:在$\mathbf{x}$处,平均权重应该反映$f$相对于$g$的``超出程度''。}

{\color{red}\textbf{[为什么这个条件保证适当加权]:} 证明思路:计算$E[h(\mathbf{X})W]$。设$(\mathbf{X}, W)$的联合密度为$g(\mathbf{x}, w)$,则
\begin{align*}
E[h(\mathbf{X})W] &= \iint h(\mathbf{x}) w \cdot g(\mathbf{x}, w) d\mathbf{x} dw \\
&= \int h(\mathbf{x}) \left[\int w \cdot g(\mathbf{x}, w) dw\right] d\mathbf{x} \\
&= \int h(\mathbf{x}) E_g(W|\mathbf{x}) g(\mathbf{x}) d\mathbf{x} \quad \text{(这里$g(\mathbf{x})$是$\mathbf{X}$的边缘密度)}\\
&= \int h(\mathbf{x}) \cdot E_g(W) \frac{f(\mathbf{x})}{g(\mathbf{x})} \cdot g(\mathbf{x}) d\mathbf{x} \quad \text{(用充要条件)}\\
&= E_g(W) \int h(\mathbf{x}) f(\mathbf{x}) d\mathbf{x} = E_g(W) \cdot E_f[h(X)].
\end{align*}
取$c = E_g(W)$即得定义!反之,若对所有$h$都有$E[h(\mathbf{X})W] = c E_f[h(X)]$,反推可得充要条件。}

{\color{red}\textbf{[经典例子]:} 标准重要抽样:$\mathbf{X} \sim g(\mathbf{x})$,权重$W = f(\mathbf{X})/g(\mathbf{X})$是$\mathbf{X}$的函数(确定性)。此时$E_g(W|\mathbf{x}) = f(\mathbf{x})/g(\mathbf{x})$ (给定$\mathbf{x}$,$W$就确定了),而$E_g(W) = \int [f(\mathbf{x})/g(\mathbf{x})] g(\mathbf{x}) d\mathbf{x} = \int f(\mathbf{x}) d\mathbf{x} = 1$。所以$E_g(W|\mathbf{x}) = 1 \cdot f(\mathbf{x})/g(\mathbf{x}) = E_g(W) \cdot f(\mathbf{x})/g(\mathbf{x})$,满足充要条件!}

{\color{red}\textbf{[权重退化问题]:} 重要抽样的主要问题:当试投密度$g$与目标密度$f$差异较大时,权重$W_i = f(\mathbf{X}_i)/g(\mathbf{X}_i)$会高度不均匀。比如99\%的样本权重接近0,只有1\%的样本有显著权重。这导致``有效样本量''远小于名义样本量$N$,估计效率很低。这就是权重退化(weight degeneracy)现象。}

在重要抽样法和标准化重要抽样法的实际应用中, 好的试抽样分布很难获得, 所以权重 \(\left\{  {{W}_{i} = f\left( {\mathbf{X}}_{i}\right) /g\left( {\mathbf{X}}_{i}\right) }\right\}\) 经常会差别很大,使得抽样样本主要集中在少数几个权重最大的样本点上。为此, 可以舍弃权重太小的样本点, 重新抽样替换这样的样本点。这种方法称为舍选控制 (Rejection Control),描述如下。

{\color{red}\textbf{[舍选控制的基本思想]:} 核心策略:1) 对权重太小($W_i < c$)的样本点,以概率$W_i/c$接受,否则舍弃并重新抽取。这类似于舍选法(rejection sampling)!2) 对被接受的样本,调整其权重,使得$({\mathbf{X}_i, W_i^*})$仍是适当加权抽样。目标:减少权重的变异性,提高有效样本量。}

首先,需要选定权重的一个阈值 \(c\) ,然后对每个样本点 \({\mathbf{X}}_{i}\) ,若 \({W}_{i} \geq  c\) ,则接受 \({\mathbf{X}}_{i}\) ; 若 \({W}_{i} < c\) ,则以概率 \({W}_{i}/c\) 接受 \({\mathbf{X}}_{i}\) ,否则舍弃 \({\mathbf{X}}_{i}\) ,重新抽取。最后,所有样本点的权重调整为新的

{\color{red}\textbf{[算法流程详解]:} 对每个从$g$抽取的$\mathbf{X}_i$,计算$W_i = f(\mathbf{X}_i)/g(\mathbf{X}_i)$:
\begin{itemize}
\item 若$W_i \geq c$:必然接受。这些是``重要''样本,目标密度$f$在此处相对较大。
\item 若$W_i < c$:以概率$W_i/c < 1$接受。权重越小,接受概率越低。若不接受,舍弃此样本,重新从$g$抽取新样本,重复此流程。
\end{itemize}
这样,权重小的样本被``稀疏化'',权重大的样本被保留,从而减少权重的不平衡。}

\[
{W}_{i}^{ * } = {p}_{c}\frac{{W}_{i}}{\min \left\{  {1,{W}_{i}/c}\right\}  } = {p}_{c}\max \left\{  {{W}_{i},c}\right\}  , \tag{3.22}
\]

{\color{red}\textbf{[公式(3.22)解释]:} 调整后的权重$W_i^* = p_c \max\{W_i, c\}$的含义:
\begin{itemize}
\item 若$W_i \geq c$:$W_i^* = p_c W_i$(乘以归一化常数)。
\item 若$W_i < c$:$W_i^* = p_c c$。注意这个样本是以概率$W_i/c$被接受的,所以其权重需要``放大''为$c$(再乘$p_c$)来补偿接受概率的降低。
\end{itemize}
直觉:被接受的低权重样本代表了多个被舍弃的样本,所以权重需要上调。等价变换:$W_i / \min\{1, W_i/c\} = W_i / (W_i/c) = c$ (当$W_i < c$时)。}

{\color{red}\textbf{[为什么需要归一化常数$p_c$]:} 由于舍选过程改变了$\mathbf{X}_i$的分布(从$g$变为$g^*$),需要$p_c$来保证$({\mathbf{X}_i, W_i^*})$仍是适当加权抽样。$p_c$本质上是接受概率的期望,即从$g$抽取一个样本被接受的概率。}

其中

\[
{p}_{c} = \int \min \left\{  {1,\frac{f\left( \mathbf{x}\right) }{{cg}\left( \mathbf{x}\right) }}\right\}  g\left( \mathbf{x}\right) d\mathbf{x} = \int \min \left\{  {g\left( \mathbf{x}\right) ,\frac{f\left( \mathbf{x}\right) }{c}}\right\}  d\mathbf{x} \tag{3.23}
\]

{\color{red}\textbf{[公式(3.23)推导]:} $p_c$的计算:从$g$抽取$\mathbf{x}$,接受概率为
\[a(\mathbf{x}) = \begin{cases} 1, & \text{if } W(\mathbf{x}) = f(\mathbf{x})/g(\mathbf{x}) \geq c \\ W(\mathbf{x})/c = f(\mathbf{x})/(cg(\mathbf{x})), & \text{if } W(\mathbf{x}) < c \end{cases} = \min\{1, f(\mathbf{x})/(cg(\mathbf{x}))\}.\]
所以总接受概率$p_c = E_{g}[a(\mathbf{X})] = \int a(\mathbf{x}) g(\mathbf{x}) d\mathbf{x} = \int \min\{1, f(\mathbf{x})/(cg(\mathbf{x}))\} g(\mathbf{x}) d\mathbf{x} = \int \min\{g(\mathbf{x}), f(\mathbf{x})/c\} d\mathbf{x}$。}

是归一化常数,如果使用标准化重要抽样法, \({p}_{c}\) 可以省略,否则, \({p}_{c}\) 可以估计为

\[
{\widehat{p}}_{c} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}\min \left\{  {1,\frac{{W}_{i}}{c}}\right\}  . \tag{3.24}
\]

这样得到的 \(\left\{  {\left( {{\mathbf{X}}_{i},{W}_{i}^{ * }}\right) ,i = 1,\ldots ,N}\right\}\) 是关于 \(f\left( \cdot \right)\) 适当加权的,被接受的 \({\mathbf{X}}_{i}\) 的分布密度为

\[
{g}^{ * }\left( \mathbf{x}\right)  = \frac{1}{{p}_{c}}\min \left\{  {g\left( \mathbf{x}\right) ,\frac{f\left( \mathbf{x}\right) }{c}}\right\}  . \tag{3.25}
\]

{\color{red}\textbf{[公式(3.25)推导]:} 接受密度$g^*(\mathbf{x})$的推导:从$g(\mathbf{x})$抽取$\mathbf{x}$,以概率$a(\mathbf{x}) = \min\{1, f(\mathbf{x})/(cg(\mathbf{x}))\}$接受。根据条件概率,被接受样本的密度为
\[g^*(\mathbf{x}) = \frac{g(\mathbf{x}) \cdot a(\mathbf{x})}{p_c} = \frac{g(\mathbf{x}) \cdot \min\{1, f(\mathbf{x})/(cg(\mathbf{x}))\}}{p_c} = \frac{\min\{g(\mathbf{x}), f(\mathbf{x})/c\}}{p_c}.\]
注意$g^*$自动满足$\int g^*(\mathbf{x}) d\mathbf{x} = 1$(因为$p_c$就是归一化常数)。}

{\color{red}\textbf{[验证适当加权性]:} 需要验证$({\mathbf{X}_i, W_i^*})$ (其中$\mathbf{X}_i \sim g^*$,$W_i^* = p_c \max\{W_i, c\}$)是适当加权抽样。计算:当$\mathbf{X} \sim g^*$时,
\begin{align*}
E[h(\mathbf{X})W^*] &= \int h(\mathbf{x}) \cdot p_c \max\{f(\mathbf{x})/g(\mathbf{x}), c\} \cdot g^*(\mathbf{x}) d\mathbf{x} \\
&= \int h(\mathbf{x}) \max\{f(\mathbf{x})/g(\mathbf{x}), c\} \cdot \min\{g(\mathbf{x}), f(\mathbf{x})/c\} d\mathbf{x}.
\end{align*}
分两种情况:1) 若$f(\mathbf{x})/g(\mathbf{x}) \geq c$,则$\max = f/g$,$\min = g$,乘积$= f$。2) 若$f(\mathbf{x})/g(\mathbf{x}) < c$,则$\max = c$,$\min = f/c$,乘积$= f$。因此$E[h(\mathbf{X})W^*] = \int h(\mathbf{x}) f(\mathbf{x}) d\mathbf{x} = E_f[h(X)]$,满足适当加权(常数$c=1$)!}

{\color{red}\textbf{[方法总结]:} 舍选控制通过``舍弃+权重调整''实现了两个目标:1) 减少权重的变异性(所有调整后权重$\geq p_c c$),2) 保持适当加权性质。代价是需要重新抽取被舍弃的样本,增加了计算量。阈值$c$的选择是关键:太小则舍弃太少,效果不明显;太大则舍弃太多,计算量大增。实践中常取权重的中位数或某个分位数。}

阈值 \(c\) 在实际中可以从权重 \(\left\{  {W}_{i}\right\}\) 选取,比如,取为 \(\left\{  {W}_{i}\right\}\) 的某个分位数。

例 3.2.3. 用 \(\mathrm{{MC}}\) 积分法计算 \(I = {\int }_{0}^{1}{e}^{x}{dx} = e - 1 \approx  {1.718}\) 。对被积函数 \(h\left( x\right)  = {e}^{x}\) 做泰勒展开得

\[
{e}^{x} = 1 + x + \frac{{x}^{2}}{2!} + \cdots
\]

取

\[
g\left( x\right)  = c\left( {1 + x}\right)  = \frac{2}{3}\left( {1 + x}\right)
\]

要产生 \(g\left( x\right)\) 的随机数可以用逆变换法,密度 \(g\left( x\right)\) 的分布函数 \(G\left( x\right)\) 的反函数为

\[
{G}^{-1}\left( y\right)  = \sqrt{1 + {3y}} - 1,0 < y < 1
\]

因此,取 \({U}_{i}\) iid \(\mathrm{U}\left( {0,1}\right)\) ,令 \({X}_{i} = \sqrt{1 + 3{U}_{i}} - 1,i = 1,2,\ldots ,N\) ,则重要抽样法的积分公式为

\[
{\widehat{I}}_{3} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}\frac{{e}^{{X}_{i}}}{\frac{2}{3}\left( {1 + {X}_{i}}\right) }
\]

渐近方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{3}\right)  = \frac{1}{N}\left( {\frac{3}{2}{\int }_{0}^{1}\frac{{e}^{2x}}{1 + x}{dx} - {I}^{2}}\right)  \approx  {0.02691}/N.
\]

如果用平均值法, 估计公式为

\[
{\widehat{I}}_{2} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{e}^{{U}_{i}}
\]

渐近方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{2}\right)  = \frac{1}{N}\left( {{\int }_{0}^{1}{e}^{2x}{dx} - {I}^{2}}\right)  \approx  {0.2420}/N \approx  {9.0} \times  \operatorname{Var}\left( {\widehat{I}}_{3}\right)
\]

是重要抽样法方差的 9 倍。

如果用随机投点法, \(h\left( x\right)  = {e}^{x} \leq  e\left( {0 < x < 1}\right)\) ,取上界 \(M = e\) ,向 \(\left\lbrack  {0,1}\right\rbrack   \times  \left\lbrack  {0,M}\right\rbrack\) 随机投点,落到 \(f\left( x\right)\) 下方的概率为

\[
p = I/\left( {M\left( {b - a}\right) }\right)  = \left( {e - 1}\right) /e,
\]

设投 \(N\) 点落到 \(h\left( x\right)\) 下方的频率为 \(\widehat{p}\) ,用随机投点法估计 \(I\) 的公式为

\[
{\widehat{I}}_{1} = \widehat{p} \cdot  M\left( {b - a}\right)  = e\widehat{p},
\]

渐近方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{1}\right)  = {e}^{2}p\left( {1 - p}\right) /N = \left( {e - 1}\right) /N \approx  {1.718}/N \approx  {7.1} \times  \operatorname{Var}\left( {\widehat{I}}_{2}\right)  \approx  {64.8} \times  \operatorname{Var}\left( {\widehat{I}}_{3}\right)
\]

可见选择合适的抽样算法对减少计算量、提高精度是十分重要的。

例 3.2.4. 设二元函数 \(f\left( {x,y}\right)\) 定义如下

\[
f\left( {x,y}\right)  = \exp \left\{  {-{45}{\left( x + {0.4}\right) }^{2} - {60}{\left( y - {0.5}\right) }^{2}}\right\}   + {0.5}\exp \left\{  {-{90}{\left( x - {0.5}\right) }^{2} - {45}{\left( y + {0.1}\right) }^{4}}\right\}
\]

求如下二重定积分

\[
I = {\int }_{-1}^{1}{\int }_{-1}^{1}f\left( {x,y}\right) {dxdy}
\]

\(f\left( {x,y}\right)\) 有两个分别以(-0.4,0.5)和(0.5, - 0.1)为中心的峰,对积分有贡献的区域主要集中在 (-0.4,0.5)和(0.5, - 0.1)附近,在其他地方函数值很小,对积分贡献很小。

用平均值法 (3.12),取点数 \(N = {10000}\) 时 \({\widehat{I}}_{2}\) 的一个估计值为 \({\widehat{I}}_{2} = {0.132}\) ,从这一次模拟估计的 \({\widehat{I}}_{2}\) 的渐近方差为 \({1.86} \times  {10}^{-5}\) 。重复模拟 \(B = {100}\) 次,得到这 100 个 \({\widehat{I}}_{2}\) 的平均值为 0.1256,样本方差为 \({1.95} \times  {10}^{-5}\) 。

用重要抽样法。取试投密度为

\[
g\left( {x,y}\right)  \propto  \widetilde{g}\left( {x,y}\right)
\]

\[
= \exp \left\{  {-{45}{\left( x + {0.4}\right) }^{2} - {60}{\left( y - {0.5}\right) }^{2}}\right\}   + {0.5}\exp \left\{  {-{90}{\left( x - {0.5}\right) }^{2} - {10}{\left( y + {0.1}\right) }^{2}}\right\}  ,
\]

\[
- \infty  < x < \infty , - \infty  < y < \infty ,
\]

这样抽取到 \(\left\lbrack  {-1,1}\right\rbrack   \times  \left\lbrack  {-1,1}\right\rbrack\) 范围外的点对积分没有贡献,因为构成 \(g\left( {x,y}\right)\) 的两个密度都很集中,所以效率损失不大。需要求使得 \(\widetilde{g}\left( {x,y}\right)\) 化为密度的比例常数。记 \(\mathrm{N}\left( {\mu ,{\sigma }^{2}}\right)\) 的分布密度为 \(f\left( {x;\mu ,{\sigma }^{2}}\right)\) ,对 \(\widetilde{g}\left( {x,y}\right)\) 积分得

\[
{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\widetilde{g}\left( {x,y}\right)  = \sqrt{{2\pi }/{90}}{\int }_{-\infty }^{\infty }f\left( {x; - {0.4},{90}^{-1}}\right) {dx} \cdot  \sqrt{{2\pi }/{120}}{\int }_{-\infty }^{\infty }f\left( {y;{0.5},{120}^{-1}}\right) {dy}
\]

\[
+ {0.5}\sqrt{{2\pi }/{180}}{\int }_{-\infty }^{\infty }f\left( {x;{0.5},{180}^{-1}}\right) {dx} \cdot  \sqrt{{2\pi }/{20}}{\int }_{-\infty }^{\infty }f\left( {y; - {0.1},{20}^{-1}}\right) {dy}
\]

\[
= \sqrt{{2\pi }/{90}}\sqrt{{2\pi }/{120}} + {0.5}\sqrt{{2\pi }/{180}}\sqrt{{2\pi }/{20}}
\]

\[
\approx  {0.1128199}
\]

于是令

\[
g\left( {x,y}\right)  = \widetilde{g}\left( {x,y}\right) /{0.1128199}
\]

\[
= {0.5358984f}\left( {x; - {0.4},{90}^{-1}}\right) f\left( {y;{0.5},{120}^{-1}}\right)
\]

\[
+ {0.4641016f}\left( {x;{0.5},{180}^{-1}}\right) f\left( {y; - {0.1},{20}^{-1}}\right) \text{ , }
\]

\[
- \infty  < x < \infty , - \infty  < y < \infty ,
\]

用复合抽样法对 \(g\left( {x,y}\right)\) 抽样, \(N = {10000}\) 时一次模拟得到的 \({\widehat{I}}_{3} = {0.126}\) ,从这一次模拟估计的 \({\widehat{I}}_{3}\) 的渐近方差为 \({6.84} \times  {10}^{-8}\) ,重复模拟 \(B = {100}\) 次,则 100 个 \({\widehat{I}}_{3}\) 的平均值为 0.1258,样本方差为 \({5.37} \times  {10}^{-8}\) 。 \({\widehat{I}}_{2}\) 的样本方差是 \({\widehat{I}}_{3}\) 的样本方差的 363 倍,如果要达到相同的估计方差, 两种方法的样本量相差三百多倍。

例 3.2.5. 标准化的重要抽样法在贝叶斯统计推断中有重要作用。例如,设独立的观测样本 \({Y}_{j}\) 服从如下的贝塔一二项分布:

\[
f\left( {{y}_{j} \mid  K,\eta }\right)  = P\left( {{Y}_{j} = {y}_{j}}\right)  = \left( \begin{matrix} {n}_{j} \\  {y}_{j} \end{matrix}\right) \frac{B\left( {{K\eta } + {y}_{j},K\left( {1 - \eta }\right)  + {n}_{j} - {y}_{j}}\right) }{B\left( {{K\eta },K\left( {1 - \eta }\right) }\right) },{y}_{j} = 0,1,\ldots ,{n}_{j}, \tag{3.26}
\]

其中 \(B\left( {\cdot , \cdot  }\right)\) 是贝塔函数, \({n}_{j}\) 为已知的正整数, \(K > 0,0 < \eta  < 1\) 为未知参数。贝塔一二项分布用于描述比二项分布更为分散的随机变量分布。按照贝叶斯统计的做法,假设参数 \(\left( {K,\eta }\right)\) 也是随机变量,具有所谓的 “先验分布”,假设 \(\left( {K,\eta }\right)\) 有如下的 “无信息” 先验分布密度:

\[
\pi \left( {K,\eta }\right)  \propto  \frac{1}{{\left( 1 + K\right) }^{2}}\frac{1}{\eta \left( {1 - \eta }\right) }, \tag{3.27}
\]

则 \(\left( {K,\eta }\right)\) 有如下的 “后验密度”:

\[
\widetilde{p}\left( {K,\eta  \mid  \mathbf{Y}}\right)  \propto  \pi \left( {K,\eta }\right) \mathop{\prod }\limits_{{j = 1}}^{n}f\left( {{y}_{j} \mid  K,\eta }\right)
\]

\[
\propto  \frac{1}{{\left( 1 + K\right) }^{2}}\frac{1}{\eta \left( {1 - \eta }\right) }\mathop{\prod }\limits_{{j = 1}}^{n}\frac{B\left( {{K\eta } + {y}_{j},K\left( {1 - \eta }\right)  + {n}_{j} - {y}_{j}}\right) }{B\left( {{K\eta },K\left( {1 - \eta }\right) }\right) }. \tag{3.28}
\]

设要求 \(E\left( {\log K \mid  \mathbf{Y}}\right)  = {\int }_{0}^{\infty }\log K\widetilde{p}\left( {K,\eta  \mid  \mathbf{Y}}\right) {dK}\) 的值。

如果可以从后验密度 \(\widetilde{p}\left( {K,\eta  \mid  \mathbf{Y}}\right)\) 直接抽样,可以用平均值法估计 \(E\left( {\log K \mid  \mathbf{Y}}\right)\) ,但从 (3.28) 来看很难直接抽样。为此,使用标准化的重要抽样法。为了解除 \(\left( {K,\eta }\right)\) 的取值限制,作变换 \(\alpha  = \log K,\beta  = \log \frac{\eta }{1 - \eta }\) ,则 \(\alpha ,\beta  \in  \left( {-\infty ,\infty }\right)\) ,而(3.28)对应的 \(\left( {\alpha ,\beta }\right)\) 的后验密度为:

\[
p\left( {\alpha ,\beta  \mid  \mathbf{Y}}\right)  \propto  \frac{{e}^{\alpha }}{{\left( 1 + {e}^{\alpha }\right) }^{2}}\mathop{\prod }\limits_{{j = 1}}^{n}\frac{B\left( {\frac{{e}^{\alpha }}{1 + {e}^{-\beta }} + {y}_{j},\frac{{e}^{\alpha }}{1 + {e}^{\beta }} + {n}_{j} - {y}_{j}}\right) }{B\left( {\frac{{e}^{\alpha }}{1 + {e}^{-\beta }},\frac{{e}^{\alpha }}{1 + {e}^{\beta }}}\right) }. \tag{3.29}
\]

取值无限制的随机变量试抽样密度经常使用自由度较小的 \(\mathrm{t}\) 分布,比如 \(\mathrm{t}\left( 4\right)\) 分布,设 \(\mathrm{t}\left( 4\right)\) 分布密度函数为 \(g\left( \cdot \right)\) ,用独立的 \(\mathrm{t}\left( 4\right)\) 分布生成 \(\left( {\alpha ,\beta }\right)\) 的试抽样样本 \(\left( {{\alpha }_{i},{\beta }_{i}}\right) ,i = 1,2,\ldots ,N\) ,可以估计 \(E\left( {\log K \mid  \mathbf{Y}}\right)\) 为

\[
\widehat{\alpha } = \frac{\mathop{\sum }\limits_{{i = 1}}^{N}{\alpha }_{i}\frac{p\left( {{\alpha }_{i},{\beta }_{i} \mid  \mathbf{Y}}\right) }{g\left( {\alpha }_{i}\right) g\left( {\beta }_{i}\right) }}{\mathop{\sum }\limits_{{i = 1}}^{N}\frac{p\left( {{\alpha }_{i},{\beta }_{i} \mid  \mathbf{Y}}\right) }{g\left( {\alpha }_{i}\right) g\left( {\beta }_{i}\right) }}. \tag{3.30}
\]

其中的 \(p\left( {{\alpha }_{i},{\beta }_{i} \mid  \mathbf{Y}}\right)\) 只要用(3.29)的右侧计算,因为分子和分母的归一化常数可以消掉。

\subsection{分层抽样法}

用平均值法计算 \({\int }_{C}h\left( \mathbf{x}\right) d\mathbf{x}\) ,若 \(h\left( \mathbf{x}\right)\) 在 \(C\) 内取值变化范围大则估计方差较大。重要抽样法选取了与 \(f\left( x\right)\) 形状相似但是容易抽样的密度 \(g\left( x\right)\) 作为试投密度,大大提高了精度,但是这样的 \(g\left( \mathbf{x}\right)\) 有时难以找到。

{\color{red}\textbf{[核心思想——分而治之]:} 分层抽样的根本思想是``divide and conquer''——将一个高方差的复杂问题分解为多个低方差的简单子问题。想象你要估计一个国家的平均收入:如果直接随机抽样,可能碰巧抽到很多富人或很多穷人,方差巨大;但如果先按地区(城市/农村)或职业分层,每层内收入相近,估计就稳定得多。数学上,这利用了方差分解:总方差 = 层内方差 + 层间方差,分层后只需估计层内小方差!}

{\color{red}\textbf{[生动比喻]:} 分层抽样就像``质量检验''——一批产品中既有合格品又有次品,如果混在一起随机检验,结果波动很大;聪明的做法是先用某种特征(如重量、外观)粗略分类,再对每类分别检验。又像``民意调查''——不同年龄段、地区的人观点差异大,直接平均会被个别极端样本主导;分层后各组内观点相近,估计更准确。}

{\color{red}\textbf{[与重要抽样对比]:} 重要抽样需要找到一个好的试投密度$g(x)$,这在高维或复杂函数时很困难——你必须对$f$的形状有精确了解。分层抽样则简单得多:只需将积分区域切分成若干块,使每块内$h(x)$变化不大即可,无需精确知道$h$的全局结构。这是``粗粒化''vs``精细模拟''的权衡。}

如果把 \(C\) 上的积分分解为若干个子集上的积分,使得 \(h\left( \mathbf{x}\right)\) 在每个子集上变化不大,分别计算各个子集上的积分再求和,可以提高估计精度。这种方法与 \(§{2.2.5}\) 的复合抽样法类似, 叫做分层抽样法。这也是抽样调查中的重要技术。

例 3.2.6. 对函数

\[
h\left( x\right)  = \left\{  \begin{array}{ll} 1 + \frac{x}{10}, & 0 \leq  x \leq  {0.5} \\   - 1 + \frac{x}{10}, & {0.5} < x \leq  1 \end{array}\right.
\]

求定积分

\[
I = {\int }_{0}^{1}h\left( x\right) {dx}
\]

{\color{red}\textbf{[精心设计的例子]:} 这个例子极具教学价值!函数$h(x)$在前半段接近+1(范围1.0到1.05),在后半段接近-1(范围-1.0到-0.95),两段``几乎相互抵消'',真实积分只有0.05——这是一个小量!如果直接平均,会采样到很多+1附近和-1附近的值,它们求和时会产生巨大的随机涨落(正负相消的不确定性),导致相对误差极大。这正是分层抽样大显身手的场景!}

可以得 \(I\) 的精确值为 \(I = {0.05}\) 。我们用平均值法和分层抽样法来估计 \(I\) 并比较精度。

在 \(\left\lbrack  {0,1}\right\rbrack\) 区间随机抽取 \(N\) 点用平均值法得 \({\widehat{I}}_{2}\) ,其渐近方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{2}\right)  = \frac{\operatorname{Var}\left( {h\left( U\right) }\right) }{N} = \frac{143}{150N} \approx  \frac{0.9533}{N}.
\]

{\color{red}\textbf{[方差来源分析]:} 这个方差为什么这么大?因为$h(U)$在$[0,1]$上的取值跨度约为2(从+1跳到-1),而期望值只有0.05。直接平均时,随机点既可能落在+1区域也可能落在-1区域,每次采样都像在做``抛硬币''——有一半机会得到+1附近的值,一半机会得到-1附近的值。这种``大幅度振荡''导致$\text{Var}(h(U)) \approx 0.95$,远大于积分真值$I=0.05$!信噪比极低。}

把 \(I\) 拆分为 \(\left\lbrack  {0,{0.5}}\right\rbrack\) 和 \(\left\lbrack  {{0.5},1}\right\rbrack\) 上的积分,即

\[
I = a + b = {\int }_{0}^{0.5}h\left( x\right) {dx} + {\int }_{0.5}^{1}h\left( x\right) {dx},
\]

{\color{red}\textbf{[分层的关键]:} 注意分层点选在$x=0.5$——这正是函数$h(x)$的跳跃点!这不是巧合,而是策略:我们把``性质差异巨大''的两个区域分开处理。第一层$h(x) \approx +1$,第二层$h(x) \approx -1$,各层内部$h$几乎不变,方差自然就小了。这就像``隔离冲突''——把意见对立的两派分开统计,避免混在一起相互干扰。}

对 \(a\) 和 \(b\) 分别用平均值法,得

\[
\widehat{a} = \frac{0.5}{N/2}\mathop{\sum }\limits_{{i = 1}}^{{N/2}}h\left( {{0.5}{U}_{i}}\right)  = \frac{0.5}{N/2}\mathop{\sum }\limits_{{i = 1}}^{{N/2}}\left( {1 + {0.05}{U}_{i}}\right) ,
\]

\[
\widehat{b} = \frac{0.5}{N/2}\mathop{\sum }\limits_{{i = \left( {N/2}\right)  + 1}}^{N}h\left( {{0.5} + {0.5}{U}_{i}}\right)  = \frac{0.5}{N/2}\mathop{\sum }\limits_{{i = \left( {N/2}\right)  + 1}}^{N}\left( {-1 + {0.05} + {0.05}{U}_{i}}\right) ,
\]

\[
{\widehat{I}}_{5} = \widehat{a} + \widehat{b},
\]

则分层抽样法结果 \({\widehat{I}}_{5}\) 的渐近方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{5}\right)  = \operatorname{Var}\left( {\widehat{a} + \widehat{b}}\right)  = \operatorname{Var}\left( \widehat{a}\right)  + \operatorname{Var}\left( \widehat{b}\right)
\]

{\color{red}\textbf{[方差相加的关键]:} 为什么$\text{Var}(\hat{a} + \hat{b}) = \text{Var}(\hat{a}) + \text{Var}(\hat{b})$?因为两层使用独立的随机点!这是分层抽样的一个微妙优势:各层估计相互独立,方差简单相加。而且注意,虽然两层的积分值几乎相互抵消($a \approx 0.525, b \approx -0.475$),但它们的估计误差是统计独立的,不会像直接平均那样产生``相消放大''效应。}

\[
= {0.25}\frac{\operatorname{Var}\left( {1 + {0.05U}}\right) }{N/2} + {0.25}\frac{\operatorname{Var}\left( {-{0.95} + {0.05U}}\right) }{N/2} = \frac{1/{4800}}{N},
\]

{\color{red}\textbf{[数学细节]:} 让我们仔细看这个计算:第一层$h(x)=1+0.05U$,方差是$\text{Var}(0.05U) = 0.05^2 \cdot \text{Var}(U) = 0.0025 \times 1/12 \approx 0.0002$。第二层同理。每层用$N/2$个样本,前面的0.25是区域长度的平方$(0.5)^2$。关键是:每层内部$h$的变化只有0.05(从1.0到1.05或从-1.0到-0.95),而不是整体的2(从+1到-1)!这就是方差缩小的根本原因。}

{\color{red}\textbf{[惊人效果]:} 分层后方差从$0.9533/N$降到$0.0002/N$,减少了约4500倍!这是一个令人震撼的数字——相同样本量下,分层抽样的标准误差缩小67倍($\sqrt{4500} \approx 67$);或者说,要达到相同精度,分层抽样只需$1/4500$的样本量!原来需要450万个样本才能达到的精度,现在1000个样本就够了。这说明原来的巨大方差几乎全部来自``层间差异''($h$在两层差距约2),分层后只剩下每层内部的微小``层内涨落''($h$在每层变化仅0.05)。}

{\color{red}\textbf{[物理直觉]:} 这个例子揭示了方差的本质:方差不是来自``总变化量'',而是来自``不可预测的随机变化''。虽然$h$在$[0,1]$上总变化很大,但这种变化是``确定性的''——我们知道前半段是+1附近,后半段是-1附近。分层后,我们把这种``已知的系统性差异''用确定性的分层结构捕捉了,只让每层内``真正随机的微小涨落''影响估计。这就是``提取信号、抑制噪声''的精髓!}

分层后的估计方差远小于不分层的结果, 可以节省样本量约 4500 倍。

一般地,设积分 \(I = {\int }_{C}h\left( \mathbf{x}\right) d\mathbf{x}\) 可以分解为 \(m\) 个不交的子集 \({C}_{j}\) 上的积分,即

\[
I = {\int }_{C}h\left( \mathbf{x}\right) d\mathbf{x} = {\int }_{{C}_{1}}h\left( \mathbf{x}\right) d\mathbf{x} + {\int }_{{C}_{2}}h\left( \mathbf{x}\right) d\mathbf{x} + \cdots  + {\int }_{{C}_{m}}h\left( \mathbf{x}\right) d\mathbf{x}
\]

{\color{red}\textbf{[一般框架]:} 分层抽样的数学结构极其清晰:将一个大积分分解为多个小积分的和。关键是这个分解既是``精确的''(各层积分和等于总积分),又是``概率独立的''(各层用不同的随机点)。这种结构确保了:1)无偏性被保留;2)方差可以独立计算后相加;3)每层可以用不同的样本量优化。这是``线性可加性''带来的优雅性质。}

在 \({C}_{j}\) 投 \({n}_{j}\) 个随机点 \({X}_{ji} \sim  \mathrm{U}\left( {C}_{j}\right) ,i = 1,\ldots ,{n}_{j}\) ,则 \(I\) 的 \(m\) 个部分可以分别用平均值法估计,由此得 \(I\) 的分层估计为

\[
{\widehat{I}}_{5} = \mathop{\sum }\limits_{{j = 1}}^{m}\frac{V\left( {C}_{j}\right) }{{n}_{j}}\mathop{\sum }\limits_{{i = 1}}^{{n}_{j}}h\left( {X}_{ji}\right)
\]

{\color{red}\textbf{[估计量结构]:} 注意每层的估计量形式:$\frac{V(C_j)}{n_j}\sum_{i=1}^{n_j} h(X_{ji})$——这就是该层的``体积''乘以``函数平均值''。整体估计是各层``加权和'',权重恰好是体积。这与直接平均$\frac{1}{N}\sum h(X_i)$的区别在于:直接平均时,样本点``随机''分布在各区域;分层后,我们``主动控制''每层的样本数,这种控制带来了方差的优化空间。}

记 \({\sigma }_{j}^{2} = \operatorname{Var}\left( {h\left( {X}_{j1}\right) }\right)\) ,划分子集时应使每一子集内 \(h\left( \cdot \right)\) 变化不大,即 \({\sigma }_{j}^{2}\) 较小。这时

\[
\operatorname{Var}\left( {\widehat{I}}_{5}\right)  = \mathop{\sum }\limits_{{j = 1}}^{m}\frac{{V}^{2}\left( {C}_{j}\right) {\sigma }_{j}^{2}}{{n}_{j}}
\]

{\color{red}\textbf{[方差公式解读]:} 这个公式是分层抽样的核心!它告诉我们总方差由各层贡献相加而成,每层的贡献是$\frac{V^2(C_j)\sigma_j^2}{n_j}$。让我们解剖它:(1)$\sigma_j^2$是层内方差——$h$在该层内的涨落;(2)$V(C_j)$是层的体积——大体积层对积分贡献大,方差也大;(3)$n_j$在分母——样本越多,方差越小。这三个因素的乘积决定了该层的``误差贡献''。优化的目标就是合理分配$n_j$使总和最小!}

若 \({\sigma }_{j}^{2}\) 可估计,应取 \({n}_{j}\) 使

\[
{n}_{j} \propto  V\left( {C}_{j}\right) {\sigma }_{j}, \tag{3.31}
\]

{\color{red}\textbf{[最优分配——Neyman分配]:} 这是统计学中著名的Neyman最优分配!公式(3.31)简洁而深刻:样本数$n_j$应该正比于$V(C_j)\sigma_j$——即``体积×标准差''的乘积。直觉上,这意味着:(1)体积大的层应该多分配样本,因为它们对总积分的贡献大;(2)方差大的层也应该多分配样本,因为它们更``难估计'',需要更多信息才能稳定。这就像``资源优化配置''——把有限的计算预算投向``重要且困难''的任务。}

{\color{red}\textbf{[与比例分配对比]:} 如果简单地按体积比例分配,$n_j \propto V(C_j)$(称为``比例分配''),也是一种常用策略。但Neyman分配更优,因为它额外考虑了$\sigma_j$——对于方差特别大的层``重点投入''。极端情况:如果某层$\sigma_j=0$(函数为常数),则该层不需要任何样本!这种``按需分配''比``平均主义''更高效。}

{\color{red}\textbf{[数学推导思路]:} 这个最优性可以用Lagrange乘数法严格证明:目标是最小化$\sum \frac{V_j^2\sigma_j^2}{n_j}$,约束是$\sum n_j = N$。构造拉格朗日函数$L = \sum \frac{V_j^2\sigma_j^2}{n_j} + \lambda(\sum n_j - N)$,对$n_j$求偏导并令其为零:$-\frac{V_j^2\sigma_j^2}{n_j^2} + \lambda = 0$,得$n_j^2 \propto V_j^2\sigma_j^2$,即$n_j \propto V_j\sigma_j$。这是一个优美的``平方根法则''!}

即

\[
{n}_{j} = N\frac{V\left( {C}_{j}\right) {\sigma }_{j}}{\mathop{\sum }\limits_{{k = 1}}^{m}V\left( {C}_{k}\right) {\sigma }_{k}},j = 1,2,\ldots ,m
\]

{\color{red}\textbf{[最优性解释]:} 为什么这样分配是最优的?因为它让``每一单位样本的贡献均等化''。在最优分配下,每层的``边际方差减少率''相同:$\frac{\partial \text{Var}}{\partial n_j} = -\frac{V_j^2\sigma_j^2}{n_j^2}$在所有层相等。这就像经济学的``边际效用均等原则''——最后一元钱花在任何商品上的效用都相同,此时总效用最大。如果某层的边际效益更高,就应该向它多分配样本,直到各层边际效益相等。}

这样取的样本量 \(\left( {{n}_{1},{n}_{2},\ldots ,{n}_{m}}\right)\) 在所有满足 \({n}_{1} + {n}_{2} + \cdots  + {n}_{m} = N\) 的取法中使得渐近方差最小 (见习题1)。

在分层抽样法中, 划分了子集后, 每一子集上的积分也可用重要抽样法计算。

{\color{red}\textbf{[组合策略]:} 这里提示了一个强大的思想——不同的方差缩减技术可以``叠加''使用!先用分层抽样将问题分解为同质性强的子问题,再在每层内用重要抽样进一步优化。这就像``先粗分、再精调''的两级优化策略,充分发挥各方法的优势。实际中,复杂问题往往需要多种技术的巧妙组合。}

分层抽样法也可以用在求随机变量函数期望的问题中。设 \(X\) 为随机变量,要求 \(X\) 的函数 \(h\left( X\right)\) 的数学期望 \(\theta  = {Eh}\left( X\right)\) 。假设存在离散型随机变量 \(Y,{p}_{j} = P\left( {Y = {y}_{j}}\right) ,j = 1,2,\ldots ,m\) , 在 \(Y = {y}_{j}\) 条件下可以从 \(X\) 的条件分布抽样,则

\[
E\left\lbrack  {h\left( X\right) }\right\rbrack   = E\{ E\left\lbrack  {h\left( X\right)  \mid  Y}\right\rbrack  \}  = \mathop{\sum }\limits_{{j = 1}}^{m}E\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{j}}\right\rbrack  {p}_{j}, \tag{3.32}
\]

{\color{red}\textbf{[条件期望的重复期望定理]:} 公式(3.32)是概率论中的``重复期望定理''(law of iterated expectations):$E[E[X|Y]] = E[X]$。这是分层抽样在随机变量情形下的理论基础。思想是:先在每个``条件世界''$Y=y_j$中计算$h(X)$的条件期望,然后对这些条件期望再按$Y$的概率分布加权平均。这种``两步平均''等价于直接对$X$平均,但分步进行给了我们优化的机会!}

{\color{red}\textbf{[选择辅助变量$Y$的艺术]:} 关键是如何选择分层变量$Y$?理想的$Y$应该满足:(1)$Y$与$h(X)$高度相关——知道$Y$后,$X$的不确定性大幅减小;(2)在给定$Y$条件下容易从$X$抽样。例如,如果$X$是收入,$Y$可以是职业、教育程度等。好的$Y$能``解释''$X$的大部分变异,使条件方差$\text{Var}[h(X)|Y=y_j]$远小于边际方差$\text{Var}[h(X)]$。}

如果在 \(Y = {y}_{j}\) 条件下生成 \(X\) 的 \({N}_{j} = N{p}_{j}\) 个抽样值,设为 \({X}_{i}^{\left( j\right) },i = 1,2,\ldots ,{N}_{j}\) ,则可以用 \(\frac{1}{{N}_{j}}\mathop{\sum }\limits_{{i = 1}}^{{N}_{j}}h\left( {X}_{i}^{\left( j\right) }\right)\) 估计 \(E\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{j}}\right\rbrack\) ,估计 \(\theta\) 为

\[
\widehat{\theta } = \mathop{\sum }\limits_{{j = 1}}^{m}\frac{1}{{N}_{j}}\mathop{\sum }\limits_{{i = 1}}^{{N}_{j}}h\left( {X}_{i}^{\left( j\right) }\right) {p}_{j} = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{m}\mathop{\sum }\limits_{{i = 1}}^{{N{p}_{j}}}h\left( {X}_{i}^{\left( j\right) }\right) , \tag{3.33}
\]

这是 \(\theta\) 的无偏和强相合估计,且估计方差

\[
\operatorname{Var}\left( \widehat{\theta }\right)  = \frac{1}{{N}^{2}}\mathop{\sum }\limits_{{j = 1}}^{m}N{p}_{j}\operatorname{Var}\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{j}}\right\rbrack
\]

\[
= \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{m}\operatorname{Var}\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{j}}\right\rbrack  {p}_{j} \tag{3.34}
\]

\[
= \frac{1}{N}E\{ \operatorname{Var}\left\lbrack  {h\left( X\right)  \mid  Y}\right\rbrack  \}  \leq  \frac{1}{N}\operatorname{Var}\left\lbrack  {h\left( X\right) }\right\rbrack  , \tag{3.35}
\]

{\color{red}\textbf{[方差缩减保证]:} 不等式(3.35)是分层抽样有效性的数学保证!它告诉我们,分层抽样的方差$\frac{1}{N}E[\text{Var}[h(X)|Y]]$永远不会超过直接平均的方差$\frac{1}{N}\text{Var}[h(X)]$。``最坏情况''是两者相等,发生在$Y$与$X$完全独立时(分层毫无意义);``最好情况''是条件方差为零,即$Y$完全决定$h(X)$(分层后无随机性)。实际中,好的分层变量能带来数量级的方差缩减!}

比直接用平均值法估计 \({Eh}\left( X\right)\) 的方差小。这里用到了条件方差的性质

\[
\operatorname{Var}\left( X\right)  = E\left\lbrack  {\operatorname{Var}\left( {X \mid  Y}\right) }\right\rbrack   + \operatorname{Var}\left\lbrack  {E\left( {X \mid  Y}\right) }\right\rbrack   \geq  E\left\lbrack  {\operatorname{Var}\left( {X \mid  Y}\right) }\right\rbrack  , \tag{3.36}
\]

{\color{red}\textbf{[方差分解定理——分层抽样的灵魂]:} 公式(3.36)是概率论中的``条件方差分解定理''(law of total variance),也是理解分层抽样威力的核心!它将总方差分解为两部分:(1)$E[\text{Var}(X|Y)]$——``层内平均方差'',即给定$Y$后$X$剩余的随机性;(2)$\text{Var}[E(X|Y)]$——``层间方差'',即不同$Y$值对应的条件期望之间的差异。}

{\color{red}\textbf{[方差分解的深刻含义]:} 让我们用具体例子理解:假设$X$是学生成绩,$Y$是班级。总方差$\text{Var}(X)$包含两部分:(1)层内方差——同一班级内学生成绩的差异(个体随机性);(2)层间方差——不同班级平均成绩的差异(系统性差异)。直接平均时,两种方差都影响估计;分层后,只有层内方差影响估计,层间方差被分层结构``吸收''了!如果班级间差异很大(层间方差大),分层的收益就巨大。}

{\color{red}\textbf{[为什么分层有效]:} 直接平均时,我们用$\frac{1}{N}\sum h(X_i)$估计$E[h(X)]$,方差是$\frac{1}{N}\text{Var}[h(X)] = \frac{1}{N}\{E[\text{Var}(X|Y)] + \text{Var}[E(X|Y)]\}$——包含两部分方差。分层抽样时,我们先算每层条件期望$E[h(X)|Y=y_j]$,再用确定性权重$p_j$加权,方差只有$\frac{1}{N}E[\text{Var}(X|Y)]$——层间方差$\text{Var}[E(X|Y)]$消失了!因为层间差异现在是``确定性地''计入,而非``随机性地''采样。这就是分层的魔力!}

如果 \(Y\) 与 \(X\) 独立则 \(E\left( {X \mid  Y}\right)  = {EX},\operatorname{Var}\left\lbrack  {E\left( {X \mid  Y}\right) }\right\rbrack   = 0\) ,这时分层抽样法比平均值法没有改进。从(3.34)可以看出,如果第 \(j\) 层样本的函数 \(h\left( {X}_{i}^{\left( j\right) }\right) ,i = 1,2,\ldots ,N{p}_{j}\) 的样本方差为 \({S}_{j}^{2}\) , 则 \(\operatorname{Var}\left( \widehat{\theta }\right)\) 的一个无偏估计是

\[
\overset{⏜}{\operatorname{Var}\left( \widehat{\theta }\right) } = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{m}{S}_{j}^{2}{p}_{j} \tag{3.37}
\]

公式(3.33)取第 \(j\) 层样本数 \({N}_{j} = N{p}_{j}\) ,仅考虑了 \(Y\) 的取值分布,而未考虑 \(X \mid  Y = {y}_{j}\) 的条件分布情况。类似于(3.31)和(3.32),应该对 \(\operatorname{Var}\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{j}}\right\rbrack\) 较大的层取更多的样本。使得估计方差最小的分层样本量分配满足 \({N}_{j} \propto  {p}_{j}\sqrt{\operatorname{Var}\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{j}}\right\rbrack  }\) ,即

\[
{N}_{j} = N\frac{{p}_{j}\sqrt{\operatorname{Var}\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{j}}\right\rbrack  }}{\mathop{\sum }\limits_{{k = 1}}^{m}{p}_{k}\sqrt{\operatorname{Var}\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{k}}\right\rbrack  }}. \tag{3.38}
\]

{\color{red}\textbf{[比例分配vs最优分配]:} 这里有两种分配策略的对比:(1)\textbf{比例分配}:$N_j = Np_j$,完全按照$Y$的概率分布分配样本,简单但未必最优;(2)\textbf{最优分配}:$N_j \propto p_j\sigma_j$(其中$\sigma_j = \sqrt{\text{Var}[h(X)|Y=y_j]}$),既考虑概率权重$p_j$,又考虑条件标准差$\sigma_j$。如果所有层的条件方差相等,两种策略一致;但如果某层条件方差特别大,最优分配会向其倾斜更多样本。}

{\color{red}\textbf{[最优分配的直觉]:} 公式(3.38)说:样本数应该正比于$p_j\sigma_j$。为什么?(1)$p_j$大的层在总期望中``权重高'',估计误差对总体影响大,应该多投入;(2)$\sigma_j$大的层``难估计'',需要更多样本才能稳定。这两个因素相乘得到最优配置。极端情况:如果某层$\sigma_j=0$(条件方差为零),那层只需1个样本即可!这体现了``按需分配''的智慧。}

{\color{red}\textbf{[与积分情形对比]:} 注意积分情形的最优分配是$n_j \propto V(C_j)\sigma_j$(公式3.31),而随机变量情形是$N_j \propto p_j\sigma_j$(公式3.38)。两者结构一致!$V(C_j)$对应$p_j$(都是``重要性权重''),$\sigma_j$都是标准差。这揭示了深层的统一性:无论是几何积分还是概率期望,最优分配的原则都是``重要性×难度''。}

在 \(\operatorname{Var}\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{j}}\right\rbrack\) 未知的时候,可以预先抽取一个小的样本估计 \(\operatorname{Var}\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{j}}\right\rbrack\) ,然后按估计的最优 \({N}_{j}\) 分配各层的样本量。采用(3.38)的分层样本量后,

\[
\widehat{\theta } = \mathop{\sum }\limits_{{j = 1}}^{m}\frac{1}{{N}_{j}}\mathop{\sum }\limits_{{i = 1}}^{{N}_{j}}h\left( {X}_{i}^{\left( j\right) }\right) {p}_{j}, \tag{3.39}
\]

\[
\operatorname{Var}\left( \widehat{\theta }\right)  = \mathop{\sum }\limits_{{j = 1}}^{m}\frac{{p}_{j}^{2}\operatorname{Var}\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{j}}\right\rbrack  }{{N}_{j}}, \tag{3.40}
\]

于是 \(\operatorname{Var}\left( \widehat{\theta }\right)\) 的估计为

\[
\widehat{\operatorname{Var}\left( \widehat{\theta }\right) } = \mathop{\sum }\limits_{{j = 1}}^{m}\frac{{p}_{j}^{2}{S}_{j}^{2}}{{N}_{j}}, \tag{3.41}
\]

其中 \({S}_{j}^{2}\) 是第 \(j\) 层样本函数 \(\left\{  {h\left( {X}_{i}^{\left( j\right) }\right) ,i = 1,2,\ldots ,{N}_{j}}\right\}\) 的样本方差。

{\color{red}\textbf{[实施策略]:} 在实践中,最优分配(3.38)需要知道各层的条件方差$\text{Var}[h(X)|Y=y_j]$。如果事先不知道,常用``两阶段策略'':(1)\textbf{导航阶段}(pilot phase):用较小样本量(如总预算的5-10\%)对各层进行初步探索,估计$\sigma_j$;(2)\textbf{主阶段}(main phase):根据导航阶段的$\hat{\sigma}_j$计算最优分配,投入剩余样本。这种``先侦察、后主攻''的策略在实际中非常有效。}

{\color{red}\textbf{[本质总结]:} 分层抽样的三大原理总结——(1)\textbf{分解原理}:把异质性大的总体分解为同质性强的子层,将``大方差问题''化为``多个小方差问题'';(2)\textbf{独立估计}:各层独立估计,层内方差小导致估计容易,且各层误差相互独立;(3)\textbf{确定性整合}:层间用确定性权重$p_j$加权求和,消除了层间方差的随机性影响。这是``化整为零、各个击破、精确整合''的完美实践!}

分层抽样法的本质是把 \(X\) 的值相近的抽样分入一层,使得同层的 \(X\) 条件方差较小,从而减小估计方差。

{\color{red}\textbf{[方法论启示]:} 分层抽样体现了处理复杂问题的一般方法论:(1)\textbf{识别结构}:发现问题中的``异质性''来源——哪些因素导致变量差异大?(2)\textbf{分而治之}:按异质性因素分层,使每层内部``同质化'';(3)\textbf{优化配置}:根据各子问题的难度和重要性合理分配资源;(4)\textbf{系统综合}:用精确的数学框架整合子问题的解。这种思想在科学计算、机器学习、工程优化中无处不在!}

例 3.2.7. 设 \(U \sim  \mathrm{U}\left( {0,1}\right)\) ,要估计 \(\theta  = {Eh}\left( U\right)  = {\int }_{0}^{1}h\left( x\right) {dx}\) 。令 \(Y = \operatorname{ceil}\left( {mU}\right)\) ,即当且仅当 \(\frac{j - 1}{m} < U \leq  \frac{j}{m}\) 时 \(Y = j,j = 1,2,\ldots ,m\) ,可以按照 \(Y\) 分层抽样估计 \(\theta\) :

{\color{red}\textbf{[均匀分层策略]:} 这个例子展示了最简单的分层方案:将$[0,1]$均匀切分为$m$个等长子区间$[\frac{j-1}{m}, \frac{j}{m}]$,每层长度$1/m$,每层概率$P(Y=j)=1/m$。这种``等距分层''虽然简单,但对很多函数已经相当有效——它保证了每层的``动态范围''不超过整体的$1/m$。对于平滑函数,这种局部小区间上的近似常数假设是合理的。}

\[
\theta  = E\left\lbrack  {h\left( U\right) }\right\rbrack   = \mathop{\sum }\limits_{{j = 1}}^{m}E\left\lbrack  {h\left( U\right)  \mid  Y = j}\right\rbrack  P\left( {Y = j}\right)  \tag{3.42}
\]

\[
= \frac{1}{m}\mathop{\sum }\limits_{{j = 1}}^{m}E\left\lbrack  {h\left( U\right)  \mid  Y = j}\right\rbrack  , \tag{3.43}
\]

易见 \(Y = j\) 条件下 \(U\) 服从 \(\left( {\frac{j - 1}{m},\frac{j}{m}}\right)\) 上的均匀分布,设 \({U}_{1},{U}_{2},\ldots ,{U}_{n}\) 是 \(\mathrm{U}\left( {0,1}\right)\) 的独立抽样, 则用分层抽样法取每层 \({N}_{j} = 1\) 估计 \(\theta  = {Eh}\left( U\right)\) 为

\[
\widehat{\theta } = \frac{1}{m}\mathop{\sum }\limits_{{j = 1}}^{m}h\left( \frac{j - 1 + {U}_{j}}{m}\right) . \tag{3.44}
\]

{\color{red}\textbf{[每层一个样本的妙处]:} 注意这里$N_j=1$——每层只取一个样本!总共用$m$个样本。估计量形式优雅:$\hat{\theta} = \frac{1}{m}\sum_{j=1}^m h(\frac{j-1+U_j}{m})$。这实际上是一种``分层随机化''策略:在第$j$层的$[\frac{j-1}{m}, \frac{j}{m}]$内随机选一点。这确保了$m$个样本点``均匀覆盖''$[0,1]$——不会像纯随机那样扎堆或留空白,这种``强制分散''大幅降低方差!}

{\color{red}\textbf{[与拟蒙特卡洛联系]:} 公式(3.44)的结构很像拟蒙特卡洛方法!如果取$U_j=0.5$(每层中点),就得到确定性的``中点法则'';如果$U_j$随机,就是分层随机化。这介于纯随机(Monte Carlo)和确定性(拟Monte Carlo)之间,兼顾了随机性(便于误差分析)和确定性(样本均匀分布)的优点。当$m$很大时,这种方法的方差衰减速度可以达到$O(1/m^2)$而非$O(1/m)$,对光滑函数尤其有效!}

{\color{red}\textbf{[分层细化]:} 层数$m$的选择权衡:(1)$m$越大,层内$h$变化越小,条件方差$\text{Var}[h(U)|Y=j]$越小,理论上越好;(2)但$m$太大时,每层样本太少,估计不稳定。实践中,$m=\sqrt{N}$(总样本$N$的平方根)常是好选择。例如$N=10000$时取$m=100$,每层100个样本。这个例子的极端情形$N_j=1$适合$h$非常光滑的情况。}
